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1.  Statement  of  Problem  Studied 


For  more  than  six  years,  Boeing  Helicopters  and  MIT  have  been  working  to  develop  actuators  for 
helicopter  rotor  control  that  can  significantly  reduce  the  vibration  typically  experienced  by  helicopters. 
This  work  has  lead  to  the  development  of  the  X-frame  actuator,  which  has  demonstrated  the  bandwidth, 
force,  and  stroke  required  for  this  application.  Boeing’s  plan  was  to  build  an  actively  controlled  rotor  with 
the  X-frame  for  an  MD-900  helicopter,  and  flight-test  the  active  rotor  system  in  2002,  as  part  of  the  Smart 
Materials  and  Structures  Demonstration  Program.  The  original  goal  of  this  project  was  to  develop  control 
approaches  and  software  suitable  for  use  in  the  flight  tests  of  the  MD-900  helicopter,  and  to  assist  Boeing 
before  and  during  the  flight  test,  to  enable  them  to  use  our  advanced  control  approaches.  Our  original 
statement  of  work  called  for  us  to  complete  the  development  of  a  system  identification  methodology  to 
allow  identification  of  the  transfer  function  of  the  rotor,  including  periodic  effects;  write  Matlab  and/or 
Simulink  software  to  implement  the  system  identification  methodology  and  continuous-time  controller,  in 
support  of  Boeing’s  flight  test  program  in  the  DARPA  Smart  Structures  and  Materials  Demonstration 
Program;  participate  in  the  flight  test  program,  by  analyzing  data  as  it  becomes  available,  and  making 
recommendations  to  the  flight  test  team.  We  would  then,  based  on  the  flight  test  data,  design  a 
continuous-time  control  law  for  flight  test.  In  addition,  we  would  evaluate  the  performance  of  the 
continuous-time  controller,  as  compared  to  the  conventional  controllers  implemented,  and  support  Boeing 
Mesa  as  needed  in  their  design,  manufacture,  and  analysis  of  the  “double  X-frame  actuator,”  previously 
developed  at  MIT. 

Due  to  delays  in  the  flight  test  test  program,  our  statement  of  work  was  modified  just  prior  to  January 
2003.  The  modification  called  for  us  to  complete  the  analysis  of  the  flight  test  data  at  Langley  on  a  model, 
active-twist  rotor,  done  jointly  with  Prof  Carlos  Cesnik  at  the  University  of  Michigan;  analyze  the  new 
X-Frame  bench  data  taken  at  Boeing,  in  order  to  improve  the  design  of  flap  position  controllers;  and 
document  the  results  of  the  Langley  experiment  and  the  bench  tests  in  archival  publications.  In  addition, 
we  would  participate  in  the  Boeing  rotor  tests  if  those  tests  if  those  tests  occur  in  the  time  frame  of  the 
effort.  (In  fact,  the  rotor  tests  did  not  occur  within  the  time  frame  of  the  grant.) 

2.  Summary  of  Important  Results 

In  this  section  we  briefly  describe  the  important  results  of  our  investigation.  More  detail  on  some 
these  results  are  given  in  the  Appendices. 

Wind  Tunnel  Test  at  Langley.  While  Boeing  prepared  for  the  flight  test,  we  collaborated  with  Prof 
Carlos  Cesnik’s  research  group  at  University  of  Michigan  to  support  his  NASA-supported  program  for 
the  Active  Twist  Rotor  (ATR)  in  forward  flight  in  the  Transonic  Dynamics  Tunnel  (TDT)  at  NASA 
Langley.  The  collaboration  was  carried  out  to  test  our  system  identification  methodology  and  continuous¬ 
time  higher  harmonic  control  (HHC)  algorithms  on  an  actual  rotor,  and  to  gain  experience  for  the 
upcoming  Boeing  tests.  For  the  wind  tunnel  test,  we  implemented  the  system  identification  algorithms  in 
Simulink  to  estimate  the  transfer  function  of  the  rotor,  and  designed  continuous-time  HHC  algorithms 
incorporating  an  anti-windup  scheme.  Some  of  the  important  results  and  conclusions  on  the  test  are  given 
below. 

•  We  analyzed  the  performance  of  our  continuous-time  controllers,  and  found  that  the  vibration  level 
was  significantly  reduced  (~  20  dB),  although  we  designed  controllers  based  on  the  assumption  that 
the  rotor  system  is  Linear  Time-Invariant  (LTI). 

•  We  applied  the  system  identification  methodologies  to  investigate  periodic  effects  in  the  rotor 
systems.  The  results  of  estimating  Harmonic  Transfer  Function  (HTF)  of  the  rotor  show  that  the 
magnitude  of  Go  (corresponding  to  LTI  transfer  function)  is  much  larger  than  other  HTF’s,  at  least 
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by  20  dB.  This  result  explains  how  we  could  achieve  20  dB  of  vibration  reduction  using 
continuous-time  controllers  based  on  LTI  plant  assumption.  This  is  important,  in  that  it  implies  that 
classical  control  laws  should  work  well  for  the  controlling  the  ATR  rotor,  which  should  greatly 
simplify  the  design  of  closed-loop  controllers  for  the  ATR  rotor.  The  following  section  will 
describe  the  system  identification  methodologies  in  more  detail. 

•  We  developed  an  anti -windup  system  for  our  controller  that  will  ensure  the  integrators  in  each  loop 
do  not  wind  up,  and  incorporated  it  in  our  continuous-time  controllers  during  the  wind  tunnel  tests. 
The  results  of  closed-loop  performance  analysis  show  that  the  anti-windup  system  worked  well  as 
expected,  and  the  closed-loop  actuation  signal  was  not  saturated. 

System  Identification  Methodology.  We  improved  our  system  identification  methodologies  to 
estimate  the  transfer  functions  of  a  linear  time  periodic  (LTP)  system.  The  transfer  properties  of  a  linear 
time  periodic  (LTP)  system  can  be  described  by  harmonic  transfer  functions  (HTF)  that  give  the  input- 
output  relationship  between  the  Fourier  coefficients  of  the  input  signal  and  those  of  the  output  signal.  The 
previous  scheme  for  identifying  harmonic  transfer  functions  of  LTP  systems  used  one  long  input  signal  to 
speed  up  the  transfer  function  estimation  process,  and  employed  a  least  square  estimation  technique  to 
obtain  HTF  estimates  [3].  Since  this  least  square  estimation  problem  is  underdetermined,  an  additional 
assumption,  that  the  transfer  function  is  smooth,  is  made  in  the  ID  scheme  to  achieve  a  well-posed 
problem.  The  estimates  are  calculated  by  applying  a  quadratic  penalty  to  the  curvature  of  the  transfer 
functions.  The  advantage  of  this  approach  is  that  we  can  generate  an  input  signal  quite  easily  to  estimate 
the  RTF’s,  and  therefore  speed  up  the  system  identification  process.  However,  this  approach  requires 
much  computer  memory  and  intensive  computation  procedure,  because  we  should  take  into  account  the 
power  spectral  density  (PSD)  and  cross-spectral  density  (CSD)  matrix  of  input  and  output  spectrum  with 
all  frequencies  simultaneously. 

In  contrast  to  the  previous  approach,  our  new  approach  uses  N  chirp  signals  to  estimate  A  harmonic 
transfer  functions.  By  doing  that,  the  PSD  matrix  of  input  spectrum  becomes  full  rank,  so  there  would  be 
no  difficulty  in  inverting  the  PSD  matrix.  In  this  method,  it  is  very  important  to  take  into  consideration 
the  time  of  application  of  each  input  relative  to  the  system  period,  in  order  to  completely  characterize  the 
system  behavior,  due  to  the  time-varying  nature  of  the  system  dynamics  during  one  period.  Some  ad-hoc 
manners  should  be  adopted  in  the  new  approach,  as  in  the  previous  one,  to  generate  ‘smooth’  transfer 
functions.  The  smoothing  process  we  adopted  for  the  new  approach  is  to  add  some  small  penalty  to  the 
PSD  matrix  of  input  spectrum  before  inverting  the  matrix.  This  smoothing  process  reduces  the 
magnitudes  of  the  RTF’s  a  little  bit,  but  reduces  noise  in  the  HTF’s  quite  significantly.  It  should  be  noted 
that  the  noise  in  the  HTF’s  is  sometimes  very  significant,  so  it  is  not  easy  to  distinguish  between  the 
HTF’s  with  negligible  magnitude  but  much  noise,  and  those  with  much  magnitude  but  significant 
variations.  However,  we  can  distinguish  them  easily  using  this  ad-hoc  manner.  That  is,  HTF’s  insensitive 
to  the  small  penalty  are  real,  while  those  sensitive  to  the  penalty  are  artifact.  In  contrast  to  the  previous 
approach,  the  new  approach  doesn’t  require  much  computer  memory  or  intensive  computation  work, 
because  we  can  compute  the  transfer  functions  at  each  frequency  independently,  without  considering  the 
whole  frequency  response  simultaneously. 

Both  identification  schemes  have  been  implemented  in  MATLAB,  and  partly  coded  in  C 
programming  language  for  maximum  computational  efficiency.  These  system  ID  methods  have  been 
validated  with  analytical  results  for  a  few  well-known  LTP  systems.  The  validation  results  show  excellent 
agreement  between  the  identified  and  analytical  transfer  functions.  We  also  applied  both  system  ID 
methods  to  investigate  periodic  effects  in  the  rotor  systems  for  the  ATR  at  NASA  Langley,  and  concluded 
that  the  periodic  effects  were  negligible. 
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Intelligent  Anti-windnp  Algorithm.  We  implemented  an  anti-windup  control  system  for  the  ATR 
wind  tunnel  test  at  NASA  Langley,  and  showed  that  it  prevented  the  closed-loop  control  signal  from 
being  saturated.  We  could  easily  implement  an  anti-windup  algorithm  for  a  single  HHC,  but  it  was  not 
straightforward  to  extend  the  algorithm  to  the  multiple  harmonic  cases.  Since  we  may  need  to  implement 
multiple  harmonic  controllers  simultaneously  to  achieve  the  performance  requirement,  it  is  necessary  that 
the  maximum  allowable  control  amplitude  for  each  anti-windup  (i.e.,  the  radius  of  a  circle  for  each  HHC) 
should  be  determined  and  distributed  in  a  reasonable  way.  The  easiest  and  simplest  way  would  be  to  use 
the  maximum  control  amplitudes  for  each  HHC  so  that  they  are  summed  to  1.  For  example,  when  there 
are  three  HHC’s  in  the  loop,  some  possible  ratios  of  the  maximum  control  amplitudes  may  include 
0.2:0. 3:0. 5,  0.3:0.3:0.4,  0.2:0.4:0.4,  etc.  This  method  guarantees  that  the  actuators  are  not  saturated. 
However,  it  is  very  likely  that  the  method  will  yield  very  conservative  control  input.  That  is,  the  actual 
control  input  to  the  actuator  may  be  much  less  than  the  maximum  allowable  control  amplitude.  The 
reason  is  that  the  control  signal  from  each  HHC  may  be  out  of  phase  in  space  (e.g.,  cosine  cyclic  vs.  sine 
cyclic  signal),  and  therefore  tend  not  to  affect  each  other.  The  experimental  results  at  Langley  showed  that 
even  the  ratio  of  0.2:0. 9:0. 9  for  IP  collective,  4P  cosine  cyclic  and  4P  sine  cyclic  signal  didn’t  saturate  the 
actuators.  This  fact  motivates  us  to  develop  the  so-called  “Intelligent  anti-wind  algorithm”  for  multiple 
HHC’s. 

There  may  be  several  ways  to  implement  anti -windup  algorithms  for  multiple  HHC’s  that  will  yield 
better  performance,  compared  with  the  conservative  method  mentioned  earlier.  Our  first  intelligent  anti¬ 
windup  algorithm  sets  up  the  initial  ratio  of  the  maximum  control  amplitudes  for  each  control  loop,  and 
increases  or  decreases  the  ratios  gradually  so  that  the  control  input  can  be  maximized  without  saturating 
the  actuators.  This  algorithm  is  quite  simple  in  that  the  initial  ratio  of  the  maximum  control  amplitudes  is 
preserved  during  the  control  operation.  However,  the  algorithm  would  yield  satisfactory  closed-loop 
results  if  the  initial  ratio  is  reasonably  selected.  We  are  also  developing  methods  to  determine  the  optimal 
ratio  of  the  maximum  control  amplitudes  for  each  HHC  in  real  time. 

The  intelligent  anti-windup  scheme  was  implemented  in  Simulink,  and  simulation  results  show  that 
the  technique  is  effective  at  preventing  integrator  windup.  The  Simulink  tool  is  available  for  future 
closed-loop  tests. 

Flap  Position  Control  Loop.  One  of  the  concerns  raised  by  Boeing  is  that  the  actuators  on  different 
blades  may  have  slightly  different  scale  factors,  so  that  the  same  voltage  input  to  all  the  actuators  may 
produce  different  flap  deflections,  and  hence  yield  an  unbalanced  input  to  the  rotor.  As  requested  by 
Boeing,  we  have  worked  on  the  feasibility  study  of  improving  the  robustness  of  X-frame  actuator  systems 
using  “inner  loop”  position  controllers  to  ensure  that  all  flaps  respond  the  same  way.  We  investigated 
various  feedback  controllers  on  the  open- loop  frequency  responses  of  X-frame  actuator  installed  in  the 
bench  test  rig,  which  were  measured  at  Boeing  in  September  2002.  The  main  result  of  our  work  is  that  we 
can  achieve  roughly  50%  of  reduction  in  the  sensitivity  of  the  plant  to  uncertainty  below  40  Hz,  without 
causing  serious  performance  degradation  above  100  Hz.  If  we  want  more  reduction  in  the  plant’s 
sensitivity  in  the  low  frequency  range,  we  have  to  accept  worse  performance  in  the  high  frequency  range, 
and  vice  versa.  Also,  we  can  achieve  quite  high  robustness  of  the  plant  at  the  harmonic  frequencies  using 
several  higher  harmonic  controllers,  as  well  as  broadband  feedback  controller. 

Our  feasibility  study  shows  that  the  inner-loop  position  controller  reduces  the  sensitivity  in  the  plant 
to  uncertainty  to  some  extent.  However,  our  recommendation  is  that  the  inner-loop  position  controller 
may  be  counter-productive,  and  therefore  should  not  be  used  unless  it  is  clear  that  the  actuator  deflections 
are  significantly  different  from  blade  to  blade,  because  the  inner-loop  position  controller  reduces  the 
stability  margins  of  the  outer-loop  vibration  control. 
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3.  Papers  submitted  or  published 

We  have  submitted  two  journal  papers  and  one  conference  paper  to  the  Journal  of  the  American 

Helicopter  Society.  A  list  of  the  papers  submitted  is  given  below. 

Manuscripts  submitted,  but  not  published: 

•  Shin,  S.  J.,  Cesnik,  C.  E.  S.,  and  Hall,  S.  R.,  “Closed-loop  Control  Test  of  NASA/ARMY/MIT 
Active  Twist  Rotor  System,  Part  I:  System  Identification,”  submitted  io  Journal  of  the  American 
Helicopter  Society. 

•  Shin,  S.  J.,  Cesnik,  C.  E.  S.,  Hall,  S.  R.,  and  Song,  K.  Y.,  “Closed- loop  Control  Test  of 
NASA/ARMY/MIT  Active  Twist  Rotor  System,  Part  II:  Controller  Implementation,”  submitted 
to  Journal  of  the  American  Helicopter  Society. 

Papers  published  in  conference  proceedings: 

•  Shin,  S.  J.,  Cesnik,  C.  E.  S.,  and  Hall,  S.  R.,  “Closed-loop  Control  Test  of  NASA/ARMY/MIT 
Active  Twist  Rotor  for  Vibration  Reduction,”  Presented  at  the  American  Helicopter  Society  59* 
Annual  Forum,  Phoenix,  Arizona,  May  6-8,  2003. 

•  S.-J.  Shin  ,  C.E.S.  Cesnik  ,  and  S.R.  Hall,  “Helicopter  Vibration  Reduction  in  Forward  Flight 
Using  Blade  Integral  Twist  Control,”  American  Helicopter  Society  58*  Annual  Forum, 
Montreal,  Canada,  June  11-13,  2002. 

Technical  reports  submitted  to  ARO: 

•  Afreen  Siddiqi  and  Steven  R.  Hall,  “Identification  of  the  Harmonic  Transfer  Functions  of  a 
Helicopter  Rotor,”  Report  AMSL  #01-01,  Active  Materials  and  Structures  Laboratory,  MIT, 
Cambridge,  MA. 

We  are  also  preparing  the  following  journal  paper  that  will  be  submitted  soon  to  the  Journal  of  Guidance, 

Control,  and  Dynamics'. 

•  Song,  K.,  Siddiqi,  A.,  and  Hall,  S.  R.,  “  Development  of  System  Identification  Methodologies  for 
Linear  Time  Periodic  Systems,”  in  preparation.  Journal  of  Guidance,  Control,  and  Dynamics. 

4.  Scientific  Personnel 

1.  Dr.  Steven  R.  Hall,  Professor  of  Aeronautics  and  Astronautics. 

2.  Ms.  Afreen  Siddiqi,  formerly  a  graduate  student  in  the  Department  of  Aeronautics  and 
Astronautics.  Earned  SM  degree  while  employed  on  project. 

3.  Dr.  Kyungyeol  Song.  Department  of  Aeronautics  and  Astronautics. 

5.  Report  of  Inventions 

None 
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A  Flap  Position  Control  Loop 


We  have  worked  on  the  feasibility  study  of  improving  the  robustness  of  X- frame  actuator  systems 
using  “inner  loop”  position  control  approaches.  This  work  was  conducted  as  requested  by  Boeing,  who 
concerned  that  the  actuators  on  different  blades  may  have  slightly  different  scale  factors,  so  that  the  same 
voltage  input  to  all  the  actuators  may  produce  different  flap  deflections,  and  hence  yield  an  unbalanced 
input  to  the  rotor.  We  have  investigated  various  feedback  control  algorithms  on  the  open-loop  frequency 
responses  of  X- frame  actuator  installed  in  the  bench  test  rig,  which  were  measured  at  Boeing  in 
September  2002,  to  ensure  that  all  flaps  respond  the  same  way.  This  section  describes  the  basic  concept, 
design  procedures,  simulation  results,  and  some  comments  for  flap  position  control  loop. 

A.l  Basic  Concept  of  Flap  Position  Control  Loop 

A  block  diagram  describing  the  open-loop  actuator  system  (i.e.,  without  the  flap  position  control 
loop)  is  shown  in  Figure  (a). 


Azimuth 


(a)  (b) 

Figure  1:  Block  diagram  of  actuator  system,  (a)  Open-loop,  (b)  Closed-loop. 

Here,  y  is  the  actual  actuator  response,  r  is  the  command  actuator  input,  and  G  is  the  system  function 
(linear  or  nonlinear)  of  the  actuator  system.  We  wanty  to  follow  r  (i.e.,  y  =  r),  especially  at  the  bandwidth 
of  interest,  so  that  the  uncertainty  in  G  does  not  affect  the  actuator  response.  As  for  the  “bandwidth  of 
interest”,  we  mean  broadband  frequency  ranges  below  6P  (i.e.  between  0  and  6P),  or  at  least  several 
harmonic  frequencies  (e.g.,  IP,  2P,  3P,  ...)  if  broadband  feedback  control  is  not  possible.  We  may  also 
use  both  broadband  and  harmonic  controllers  at  the  same  time,  in  order  to  satisfy  the  requirement  of 
reducing  the  effect  of  the  uncertainty  on  the  actuator  response.  In  any  case,  the  block  diagram  of  closed- 
loop  actuator  system,  including  the  flap  position  control  loop,  can  be  represented  as  in  Figure  (b).  Using 
the  flap  position  control  loop,  the  actuator  response  y  is  expected  to  follow  its  command  input  r  within  the 
bandwidth  of  interest,  regardless  of  the  uncertainty  in  G. 

A.2  Transfer  Function  of  Actuator  Systems 

Before  designing  any  controller,  it  is  necessary  to  estimate  the  open-loop  transfer  function  of  the 
system  and  investigate  their  properties.  Figure  1  shows  the  open-loop  frequency  responses  of  the  X- frame 
actuator  installed  in  the  bench  test  rig,  which  were  measured  in  September  2002  at  Boeing.  Note  that  the 
frequency  responses  in  Figure  1  include  the  dynamics  of  high-gain  amplifier.  Also,  note  that  a  2  kHz  4- 
pole  Bessel  filter,  which  was  used  during  the  data  acquisition  probably  to  reduce  the  sensor  noise  and 
aliasing  effect,  was  replaced  by  a  4  kHz  Butterworth  filter  to  facilitate  the  design  of  feedback  controllers. 

Frequency  responses  in  Figure  1  show  some  undesirable  features  of  the  actuator  dynamics  for 
controller  design.  First,  the  actuator  response  has  a  large  resonance  around  100  Hz,  with  varying 
frequency  depending  on  the  boundary  condition,  and  several  other  resonance  above  200  Hz.  Recalling 
that  the  bandwidth  of  interest  for  our  applications  is  below  40  Hz  (~  6P),  the  effects  of  the  resonance 
around  100  Hz  on  the  controller  design  would  be  significant.  The  effect  of  resonance  could  be  reduced 
using  a  notch  filter,  if  they  do  not  vary  considerably  with  frequency.  However,  as  shown  in  Figure  1,  the 
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frequency  of  the  resonance  in  our  actuator  system  changes  quite  significantly,  which  implies  that  a  notch 
filter  cannot  be  used  to  cancel  it  out.  Another  undesirable  feature  of  the  actuator  dynamics  is  the 
substantial  phase  lag  in  the  response  above  300  Hz,  especially  in  the  transfer  function  with  light  color 
(clamp  position  6).  The  phase  lag  limits  seriously  the  magnitude  of  controller,  and  therefore  the 
performance  of  closed-loop  systems;  the  loop  gain  should  be  much  less  than  unity  above  300  Hz  to  avoid 
any  instability  problem,  which  also  limits  the  loop  gain  below  the  bandwidth  of  interest  (~40  Hz). 
Therefore,  the  actuator  system  we  are  dealing  with  forces  us  to  accept  a  trade-off  between  high 
performance  below  40  Hz  (i.e.,  high  robustness  to  uncertainty)  and  low  performance  above  300  Hz  (i.e., 
high  sensitivity  to  uncertainty). 


10“  10'  10^ 


Frequency  (FIz) 

Figure  1:  Transfer  functions  of  the  X-frame  actuator  with  a  high-gain  amplifier. 


A.3  Design  of  Controllers  and  Simulation  Results 

In  order  to  investigate  the  feasibility  of  applying  flap  position  control  loop,  we  designed  various 
feedback  controllers  and  compared  the  advantages  and  disadvantages  of  each  one.  Here,  we  present  the 
design  concept  and  simulation  results  for  three  controllers  among  them. 

(1)  The  first  controller:  classical  (non  model-based)  &  broadband  controller 

In  the  first  controller  design  approach,  we  adopted  the  classical  control  design  technique  (i.e.,  non 
model-based  approach),  aiming  at  achieving  a  ‘good  loop  shape’  in  a  broadband  frequency  range.  Here, 
‘good  loop  shape’  means  that  the  loop  transfer  function  will  have  as  high  loop  gain  as  possible  below  40 
Hz,  while  stabilizing  the  closed-loop  system  with  enough  gain/phase  margins.  It  was  our  belief  that  the 
loop  gain  should  be  at  least  1 0  dB  within  the  bandwidth  of  interest,  in  order  to  achieve  the  closed-loop 
system  with  acceptable  performance  (note  that  the  loop  gain  of  10  dB  will  reduce  the  sensitivity  of  the 
actuator  to  uncertainty  by  about  3  times,  i.e.,  ~  70%).  However,  it  turned  out  that  this  level  of  loop  gain 
could  not  be  achieved  without  serious  stability  problems  in  the  high  frequency  ranges,  due  to  the 
undesirable  dynamics  described  in  the  previous  section.  After  some  trade-off  analysis  between  the 
performance  within  the  bandwidth  of  interest  (below  40  Hz)  and  stability  problems  in  the  high-frequency 
ranges  (above  1 00  Hz),  we  concluded  that  the  maximum  loop  gain  we  could  achieve  at  40  Hz  using  a 
classical  approach,  without  compromising  the  gain/phase  margins,  is  about  5  dB. 

The  frequency  response  of  the  first  controller  is  shown  in  Figure  2(a),  and  the  resulting  loop  transfer 
function  and  sensitivity  transfer  function  for  each  actuator  response  (clamp  position  1  and  6)  are  shown  in 
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Figure  2(b).  Here,  the  first  and  second  plant  corresponds  to  the  transfer  function  with  dark  (clamp  position 
1)  and  light  (clamp  position  6)  color  in  Figure  1,  respectively.  Also,  note  that  the  sensitivity  transfer 
function  provides  the  information  about  the  reduction  of  the  closed-loop  actuator  system’s  sensitivity  to 
uncertainty,  compared  with  its  open-loop  counterpart.  In  other  words,  if  the  magnitude  of  the  sensitivity 
transfer  function  is  -20  dB,  it  means  that  the  sensitivity  of  the  actuator  response  to  uncertainty  is  reduced 
by  20  dB  (i.e.,  10  times)  using  the  feedback  controller.  The  loop  transfer  functions  in  Figure  2(b)  show 
that  the  first  controller  actually  achieves  about  6  dB  of  loop  gain  below  40  Hz,  with  acceptable  gain  and 
phase  margins.  Therefore,  we  can  expect  that  the  sensitivity  of  the  closed-loop  actuator  system  to 
uncertainty  will  be  reduced  by  a  factor  of  2  below  40  Hz.  The  expectations  are  supported  by  the 
sensitivity  transfer  functions  shown  in  Figure  2(b),  whose  magnitude  is  lower  than  -6  dB  below  40  Hz. 

As  mentioned  earlier,  this  level  of  performance  was  a  result  of  many  trade-off  analyses  between  the 
performance  and  stability  in  the  low-  and  high-frequency  ranges,  respectively.  If  we  want  to  achieve 
better  performance  (higher  loop  gain)  below  40  Hz,  we  have  to  accept  worse  stability  (lower  gain/phase 
margin)  above  1 00  Hz. 

The  closed-loop  transfer  function  for  each  actuator  response  is  shown  in  Figure  3  and  Figure  4, 
respectively,  with  both  zoomed-out  (a)  and  zoomed-in  (b)  versions.  In  order  to  see  how  robust  the  closed- 
loop  actuator  system  is  to  uncertainty,  the  magnitude  of  the  open-loop  actuator  response  was  increased 
and  decreased  by  30  %  artificially,  and  the  closed-loop  performance  of  each  case  was  compared  with  that 
of  the  nominal  one.  The  results  of  Figure  3  and  Figure  4  indicate  that  the  variation  in  the  magnitude  of  the 
closed-loop  actuator  system  with  uncertainty  is  within  around  1.5  dB  below  40  Hz.  Considering  that  30  % 
of  uncertainty  corresponds  to  roughly  3  dB,  we  can  conclude  that  the  first  controller  reduced  the 
sensitivity  of  the  actuator  at  least  by  a  factor  of  2  below  40  Hz. 

(2)  The  second  controller:  classical  &  broadband/harmonic  controller 

The  first  controller  aimed  at  improving  the  robustness  of  the  actuator  in  a  broadband  frequency  range, 
i.e.  between  0  and  40  Hz.  As  described  earlier,  there  was  a  limit  in  increasing  the  gain  for  the  broadband 
feedback  controller  (and  therefore,  improving  the  performance  of  the  closed-loop  system),  due  to  the 
substantial  resonant  peak  and  phase  delay  in  the  actuator  response.  Considering  these  limitations,  we 
designed  the  second  controller,  which  wraps  several  higher  harmonic  control  loops  around  the  closed- 
loop  actuator  system  with  the  first  controller.  In  other  words,  the  second  controller  contains  the  first 
controller  and  several  higher  harmonic  controllers.  Using  the  proposed  controller,  the  actuator  response 
follows  its  command  input  at  each  harmonic  frequencies  1P~6P  regardless  of  the  uncertainty  in  the 
system. 

The  frequency  responses  of  the  controller,  the  resulting  loop  transfer  functions  and  sensitivity  transfer 
functions,  and  the  closed-loop  actuator  transfer  functions  are  shown  in  Figure  5  -  Figure  7.  Basically,  all 
of  the  results  shown  here  are  the  same  as  those  obtained  using  the  first  controller,  except  at  each  harmonic 
frequencies  1P~6P  (here,  IP  is  assumed  to  be  6.5  Hz).  That  is,  the  second  controller  reduces  the 
sensitivity  of  the  actuator  to  uncertainty  by  a  factor  of  2  below  40  Hz  except  the  harmonic  frequencies, 
without  making  the  performance  above  1 00  Hz  worse  too  much.  At  the  harmonic  frequencies,  the  loop 
gains  are  infinite,  and  therefore  the  magnitude  of  the  sensitivity  transfer  function  becomes  zero,  which 
implies  that  the  uncertainty  in  the  system  or  environment  does  not  affect  the  actuator  response  at  the 
harmonic  frequencies  at  all.  These  phenomena  can  be  observed  from  the  closed-loop  responses  in  Figure 
6  and  Figure  7. 
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Figure  2:  Frequency  response  of  the  first  controller  (a)  and  its  loop  transfer  function  (b). 
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Figure  3:  Closed-loop  transfer  function  of  the  first  plant  obtained  using  the  first  controller. 
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Figure  4:  Closed-loop  transfer  function  of  the  second  plant  obtained  using  the  first  controller. 
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Figure  5:  Frequency  response  of  the  second  controller  (a)  and  its  loop  transfer  function  (b). 
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Figure  6:  Closed-loop  transfer  function  of  the  first  plant  obtained  using  the  second  controller. 
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Figure  7:  Closed-loop  transfer  function  of  the  second  plant  obtained  using  the  second  controller. 
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(3)  The  third  controller:  modern  (model-based)  &  broadband  controller 

The  first  and  second  controllers  were  designed  using  the  classical  control  design  approach,  where 
controllers  are  designed  based  on  the  measured  transfer  function  without  using  any  model.  In  the  third 
controller  design,  we  adopted  a  ‘modem  approach’,  in  that  a  plant  model  is  identified  from  the  measured 
transfer  function,  and  a  controller  is  designed  on  the  identified  model.  We  also  employed  a  ‘robust’ 
control  design  technique,  considering  that  the  actuator  dynamics  may  change  depending  on  its 
environment  and  that  the  controlled  system  should  be  able  to  function  adequately  regardless  of  the 
environmental  uncertainty.  Here,  we  assumed  that  there  are  two  main  uncertainties  in  the  actuator  system, 
which  are  the  first  resonant  frequency  and  the  DC  gain  (or,  the  scale  factors).  The  controller  was  designed 
using  the  fj,  synthesis  technique  (D-K  iteration),  so  that  the  controlled  actuator  system  can  tolerate  15%  of 
uncertainty  in  both  the  first  resonant  frequency  and  the  DC  gain,  without  degrading  the  closed-loop 
performance. 

Note  that  we  had  to  assume  that  the  uncertainties  in  the  resonant  frequency  and  gain  are  complex, 
although  they  are  actually  real,  because  the  current  MATLAB  version  does  not  support  real  fj,  synthesis 
(although  it  supports  real  jU  analysis).  Since  the  controller  obtained  using  the  complex  pc  synthesis  may  be 
conservative,  we  iterated  the  controller  design  via  D-K  iteration  and  used  the  real  pi  analysis  technique 
simultaneously,  to  check  the  robust  stability  and  performance  and  improve  the  closed-loop  performance. 

Figure  8  shows  the  frequency  response  of  measured  and  identified  (i.e.,  model)  transfer  function  of 
the  actuator  system.  The  frequency  responses  of  the  controller,  the  resulting  loop  transfer  functions,  and 
the  sensitivity  transfer  functions  are  shown  in  Figure  9  -  Figure  11.  The  loop  transfer  function  in  Figure 
9(b)  show  that  the  third  controller  actually  achieves  at  least  20dB  of  loop  gain  below  40  Hz,  which  is 
about  5  times  higher  than  the  first  controller,  with  acceptable  gain  and  phase  margins.  This  indicates  that 
the  sensitivity  of  the  closed-loop  actuator  system  to  uncertainty  can  be  reduced  at  least  by  a  factor  of  1 0 
below  40  Hz.  The  closed- loop  transfer  functions  for  each  actuator  response,  which  are  shown  in  Figure  10 
and  Figure  11,  respectively,  support  this  expectation.  The  results  in  Figure  10  and  Figure  1 1  show  that  the 
variation  in  the  magnitude  of  the  closed-loop  actuator  system  with  uncertainty  is  almost  negligible  below 
40  Hz.  This  is  a  significant  improvement  of  the  closed-loop  performance,  compared  with  the  first 
controller. 

However,  it  should  be  reminded  that  the  data  we  used  to  identify  the  actuator  transfer  function  and 
design  the  controller  has  frequency  components  up  to  400  Hz;  we  do  not  have  any  information  about  the 
dynamic  behavior  of  the  actuator  above  that  frequency.  Depending  on  the  high-frequency  dynamics  of  the 
actuator  (which  we  do  not  know  currently),  the  results  shown  here  may  be  changed  significantly. 


Figure  8:  Measured  and  identified  transfer  function  of  the  actuator. 


14 


Phase  (deg)  Magnitude  (dB)  ,  P^^^se  (deg)  Magnitude  (dB)  ,  Phase  (deg) 


Figure  9:  Frequency  response  of  the  third  controller  (a)  and  its  loop  transfer  function  (h). 


Norminal  plant 

■  0.7  times  nominal  plant 

■  1.3  times  nominal  plant 


10" 


10' 


_ 

Norminal  plant 

—  0.7  times  nominal  plant 

—  1.3  times  nominal  plant 

Frequency  (Hz) 


(a) 


(b) 


Figure  10:  Closed-loop  transfer  function  of  the  first  plant  obtained  using  the  third  controller, 
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Figure  11:  Closed-loop  transfer  function  of  the  second  plant  obtained  using  the  third  controller. 
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A.4  Discussions 

As  requested  by  Boeing,  we  have  worked  on  the  feasibility  study  of  making  X- frame  actuator  systems 
more  robust  using  “inner  loop”  position  controllers,  to  ensure  that  all  flaps  respond  the  same  way.  We 
investigated  various  feedback  controllers  on  the  open-loop  frequency  responses  of  X- frame  actuator 
installed  in  the  bench  test  rig,  which  were  measured  at  Boeing  in  September  2002,  in  order  to  improve  the 
robustness  of  the  actuator  system.  The  main  conclusion  of  our  work  is  that  we  can  achieve  roughly  50% 
of  reduction  in  the  sensitivity  of  the  actuator  to  uncertainty  below  40  Hz  using  the  classical  control 
approach,  and  at  least  90  %  of  reduction  using  the  modern  control  approach,  without  causing  serious 
performance  degradation  above  100  Hz.  If  we  want  more  reduction  in  the  actuator’s  sensitivity  in  the  low 
frequency  range,  we  have  to  accept  worse  performance  in  the  high  frequency  range,  and  vice  versa.  Also, 
we  can  achieve  quite  high  robustness  of  the  actuator  system  at  the  harmonic  frequencies  using  several 
higher  harmonic  controllers,  as  well  as  broadband  feedback  controllers. 

Although  the  results  of  the  feasibility  study  for  the  flap  position  control  look  promising,  especially  in 
the  third  controller  (robust&  model-based),  the  actual  performance  of  X- frame  actuator  systems  with  the 
controller  may  not  be  as  good  as  the  simulation  results,  or  even  counter-productive,  due  to  the  following 
reasons: 

(1)  As  mentioned  briefly  before,  we  do  not  have  any  information  about  the  dynamic  behavior  of  the 
actuator  above  400  Hz.  The  data  we  used  to  investigate  the  feasibility  of  the  flap  position  control, 
which  were  taken  at  Boeing  in  September  2002,  have  frequency  components  up  to  400  Hz. 
Therefore,  the  actual  performance  may  be  seriously  worse  depending  on  the  high-frequency 
dynamics  of  the  actuator,  especially  if  the  actuator  system  has  some  undesirable  dynamics  around 
400Hz,  such  as  high  resonance  and  phase  delay. 

(2)  In  order  to  achieve  20  dB  of  loop  gain  at  40  Hz,  which  will  reduce  the  sensitivity  of  the  actuator 
system  to  uncertainty  by  90%,  the  bandwidth  of  the  actuator  should  be  at  least  400  Hz,  which 
seems  significantly  higher  than  the  bandwidth  of  the  current  X- frame  actuator. 

Considering  these  facts,  our  recommendation  is  that  the  flap  position  control  may  be  counter¬ 
productive,  and  therefore  should  not  be  used  unless  it  is  clear  that  the  actuator  deflections  are  significantly 
different  from  blade  to  blade. 

B  System  Identification  of  Linear  Time  Periodic  Systems 

We  have  improved  our  system  identification  methodologies  to  estimate  the  transfer  functions  of  a 
linear  time  periodic  (LTP)  system.  In  spite  of  their  importance,  LTP  systems  have  had  much  less  attention 
than  Linear  Time  Invariant  (LTI)  systems,  primarily  due  to  the  complexity  in  analyzing  LTP  systems. 
Understanding  the  dynamic  behavior  of  LTP  systems  is  important  in  many  engineering  applications, 
especially  the  rotor  system,  in  order  to  design  effective  control  algorithms  and  improve  the  performance 
for  those  systems. 

One  of  the  fundamental  differences  between  LTP  and  LTI  systems  remains  in  the  frequency 
responses  of  those  systems  to  a  sinusoidal  input.  For  LTI  systems,  a  sinusoidal  input  at  a  given  frequency 
produces  the  system  output  with  the  same  frequency,  with  possibly  different  magnitude  and  phase.  On  the 
other  hand,  if  a  sinusoidal  signal  is  input  to  LTP  systems,  possibly  several  (or  an  infinite  number  of) 
harmonics  will  appear  in  the  output  signal,  each  with  possibly  different  magnitude  and  phase  [1].  Due  to 
these  reasons,  many  analysis  and  controller  design  tools  that  have  been  developed  for  LTI  systems,  such 
as  the  concept  of  transfer  functions.  Bode  and  Nyquist  theory,  etc.,  cannot  be  directly  applied  to  LTP 
systems. 

This  section  describes  our  system  identification  methodologies  for  LTP  systems,  which  were 
developed  based  on  the  work  of  Wereley  [2].  The  section  is  organized  as  follows.  In  Section  B.l,  an 
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overview  of  harmonic  transfer  functions  and  important  properties  of  LTP  systems  are  presented,  along 
with  a  brief  discussion  on  the  development  of  the  harmonic  state-space  model.  Section  B.2  deals  with  the 
theoretical  development  of  two  identification  schemes  we  developed  for  determining  harmonic  transfer 
functions  of  LTP  systems.  The  assumptions,  simplifications,  and  noise  reduction  techniques  of  the 
schemes  are  also  discussed.  Section  B.3  presents  the  validation  results  of  our  identification  schemes  by 
analyzing  an  LTP  system,  and  comparing  the  transfer  functions  obtained  analytically  and  empirically 
through  our  algorithms. 

B.l  Frequency  Response  of  LTP  Systems 

In  this  section,  we  review  briefly  the  characteristics  of  linear  time  periodic  (LTP)  systems  and  a  linear 
operator  to  describe  transfer  properties  of  LTP  systems.  In  LTP  systems,  the  coefficients  of  the 
differential  equations  that  describe  the  dynamics  are  time-varying  and  periodic.  These  systems  are 
described  by  a  state  space  model  of  the  form 

x(t)=A(t)x(t)  +  B(t)u(t),  y{t)  =  C(t)x(t)  +  D(t)u(t)  (1) 

where  the  matrices  A(t),  B(t),  C(t),  and  D(t)  are  periodic  with  period  T.  In  other  words, 

A(t  +  NT)  =  A(0,B(I  +  NT)  =  B(0,C(f  +  NT)  =  C(0,D(f  -i-  NT)  =  D(0  (2) 

for  any  integer  N. 

There  is  a  fundamental  difference  in  the  frequency  response  between  linear  time  invariant  (LTI) 
systems  and  LTP  systems.  For  LTI  systems,  the  response  to  an  input  sinusoid  at  a  single  frequency  is 
another  sinusoid  at  the  same  frequency  as  the  input  signal,  with  possible  changes  in  magnitude  and  phase. 

On  the  other  hand,  when  a  complex  exponential  (or  sinusoid)  e""'  is  used  to  excite  an  LTP  system,  the 
output  response  consists  of  a  superposition  of  sinusoids  not  only  at  the  input  frequency  o),  but  also  at 
several  (possibly  an  infinite  number)  other  frequencies,  co+ncjp,  each  with  its  own  magnitude  and  phase 
[1],  where  n  is  an  integer,  and  u)p=2'KlT  is  the  pumping  frequency.  Since  the  frequencies  nWp  are 
harmonics  of  the  pumping  frequency,  the  frequencies  co  +na)p  are  shifted  harmonics,  and  we  often  refer  to 
these  frequencies  simply  as  “harmonics.” 

As  expected  from  the  difference  in  the  frequency  response,  a  complex  exponential  e‘,  which  is  an 
eigenfunction  of  LTI  systems,  cannot  be  an  eigenfunction  of  LTP  systems,  because  it  results  in  a  one  to 
many  map.  In  order  to  find  eigenfunctions  of  LTP  systems,  Wereley  et  al.  defined  a  complex 
exponentially  modulated  periodic  (BMP)  function,  which  is  expressed  as  the  complex  Fourier  series  of  a 
periodic  signal  of  frequency  (Op,  modulated  by  a  complex  exponential  signal  [1], 

«GZ 

where  =  s  +  jnco^  (s  G  C) .  They  showed  that  BMP  functions  are  eigenfunctions  of  LTP  systems, 

which  implies  that  an  LTP  system  induces  at  steady  state  a  one-to-one  map  from  BMP  inputs  to  BMP 
outputs  of  the  same  frequency,  but  with  possibly  different  magnitude  and  phase  [1]. 

We  should  apply  an  BMP  function  to  LTP  systems  to  figure  out  their  transfer  properties,  as  we  do  for 
LTI  systems  using  a  complex  exponential  function.  Since  the  steady  state  response  of  an  LTP  system  in 
Bquation  1  to  an  BMP  function  is  also  an  BMP  function,  we  can  represent  the  state  and  output  variables  as 

=  ^xy"‘,y(t)  =  ^yy"‘  (4) 

«GZ  ttGZ 

Also,  the  time  periodic  coefficient  matrix  in  Bquation  1  can  be  expanded  in  a  complex  Fourier  series  as 

A(0  -  2  .C(,)  -  p) 

mG.Z  mGZ  mG.Z  mG.Z 
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Using  Equation  4  and  5,  and  applying  the  harmonic  balance  approach,  the  state  space  model  in  Equation  1 
can  be  transformed  in  the  frequency  domain  as 

(6) 

mE.Z  mE.Z  mG.Z  mE.Z 

After  expressing  the  summations  in  matrix  form,  the  equations  reduce  to 

sX  =  (A-N)X  +  BU,Y  =  CX  +  DU  (7) 

The  state  vector  X  represents  the  states  at  various  harmonics  of  a  given  frequency,  and  is  defined  as 

^^  =[•••  X-2  X-i  *0  xf  4  •••]  (8) 

The  control  signal,  U,  and  output  signal,  Y,  are  similarly  defined.  The  dynamics  matrix  A  is  a  doubly- 
infinite  Toeplitz  matrix,  given  by 


The  matrix  A„  is  the  nth  Fourier  coefficient  of  A{t).  The  matrices  B,  C,  and  D  are  similarly  defined.  The 
modulation  frequency  matrix,  N,  is  an  infinite  block  diagonal  matrix,  given  by 


0 

0-jm^I  (10) 

1- J%I 

where  I  is  the  identity  matrix  of  the  same  dimensions  as  that  of  A„. 

The  harmonic  transfer  function  (HTF)  is  defined  as  an  operator  that  relates  harmonics  of  the  input 
signal  to  harmonics  of  the  output  signals,  and  is  given  by  [1] 

G(5)  =  C[5l-(A-A^)]'‘5  +  D  (11) 

The  LTP  transfer  function  is  thus  analogous  to  the  widely  used  LTI  transfer  function,  in  that  it 
describes  the  input-output  properties  of  LTP  systems  in  the  frequency  domain.  It  should  be  noted  that  the 
expression  for  harmonic  transfer  functions  of  LTP  systems  in  Equation  1 1  has  a  very  similar  form  as  that 
for  transfer  functions  of  LTI  systems.  Although  the  matrices  in  Equation  7  are  infinite  (due  to  the  infinite 
Fourier  coefficients),  for  practical  purposes  we  can  truncate  the  number  of  terms  in  the  Fourier  series,  and 
simply  use  the  smallest  number  of  harmonics  that  adequately  represent  the  system  dynamics.  Since  the 
harmonics  generally  get  smaller  with  increasing  harmonic  number,  only  the  consideration  of  the  first  few 
harmonics  is  usually  adequate  in  describing  the  system  behavior. 

We  can  now  employ  this  notion  of  harmonic  transfer  functions  in  developing  methods  for  system 
identification  of  LTP  systems. 
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B.2  System  Identification  Schemes  for  LTP  systems 

In  this  section,  we  present  our  approaches  to  identify  harmonic  transfer  functions  of  LTP  systems.  We 
first  illustrate  the  difficulties  encountered  when  identifying  LTP  systems,  and  the  two  different 
approaches  we  developed  to  overcome  these  difficulties  will  be  described  in  detail. 

Challenges  in  LTP  System  Identification 

In  an  LTP  system,  an  input  sinusoid  at  a  single  frequency  generates  a  superposition  of  sinusoids  at 
infinite  frequencies  of  various  magnitudes  and  phase  in  the  output.  This  frequency  behavior  of  LTP 
systems  is  fundamentally  different  from  that  of  LTI  systems,  where  an  input  sinusoid  at  a  single 
frequency  generates  an  output  sinusoid  at  the  same  frequency  as  the  input  signal.  Due  to  the  different 
behaviors  in  the  frequency  domain  between  LTP  and  LTI  systems,  the  method  of  estimating  LTI  transfer 
functions,  which  is  well  established,  cannot  be  used  to  estimate  harmonic  transfer  functions  (HTF)  of  an 
LTP  system. 

To  illustrate  the  difficulty  of  identifying  LTP  systems,  we  will  initially  account  for  only  three 
frequencies  in  the  output  of  an  LTP  system,  for  each  input  frequency.  The  LTP  system  will  have  period  T, 
and  corresponding  fundamental  frequency  o)p.  We  therefore  assume  that  the  output,  T,  at  each  frequency 
CO,  comprises  of  the  linear  combination  of  the  responses  due  to  inputs  at  frequencies  co,  co  +C0p  and  co  -cOp. 

The  system  output  can  then  be  assumed  to  be  a  linear  combination  of  three  different  transfer  functions 
(each  corresponding  to  one  of  the  three  frequencies  in  the  output):  Go,  Gi,  and  G.i,  respectively.  Y  can 
thus  be  expressed  as 

Y[jo})  =  Go(jm)f/(jm)  +  G^{ja})u[ja}  -  +  joj^)  (12) 

Equivalently  in  the  time  domain,  the  output  may  be  written  as 

y{^)  =  .?o(0  *  “(0  +  ^i(0  *  [“(0^'“"']  +  ^-i(0  *  [“(0^"'“'']  (13) 

where  gk{t)  is  the  Afh  impulse  response  of  the  system,  and  denotes  the  convolution  operator. 

From  Equation  13,  one  can  see  that  the  output  y(t)  is  defined  as  the  total  response  due  to  an  input  u{t) 
that  has  been  modulated  appropriately,  and  then  convolved  with  the  respective  impulse  response 
functions.  We  can  state  more  specifically  that  the  «th  transfer  function  is  defined  to  be  the  linear  operator 
that  maps  the  output  at  frequency  co  to  an  input,  at  frequency  co,  modulated  with  This  linear  system 

can  be  represented  in  the  block  diagram  of  Figure  12. 


Figure  12:  LTP  system  model  with  three  transfer  functions. 

As  is  evident  from  Equations  12  and  13,  we  have  three  transfer  functions.  Go,  Gi,  and  G.\,  that  need  to 
be  estimated.  For  a  given  input  G(/'co)  and  resulting  output  T(/co),  we  thus  have  three  unknowns,  but  only 
one  equation.  Our  identification  problem  is  therefore  underdetermined.  In  order  to  overcome  this 
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difficulty,  we  have  developed  two  different  approaches  that  allow  us  to  estimate  RTFs  of  an  LTP  system 
systematically.  The  following  subsections  describe  our  system  identification  methods  for  LTP  systems. 


First  method 

The  basic  idea  of  the  first  method  is  to  assume  that  the  transfer  functions  are  “smooth”  and  to 
formulate  the  original  problem  as  the  minimization  of  a  cost  function  that  penalizes  the  quadratic 
estimation  error  and  the  curvature  of  harmonic  transfer  functions  (RTFs)  [3].  Rere,  “smooth”  transfer 
functions  mean  that  there  are  no  rapid  variations  with  frequency  in  the  response.  By  using  this  approach, 
the  originally  underdetermined  problem  becomes  determinate,  and  therefore,  the  RTFs  of  LTP  systems 
can  be  estimated  without  any  difficulty.  In  order  to  estimate  the  transfer  functions  with  enough  accuracy, 
the  input  signal  used  to  excite  the  system  should  have  enough  frequency  components,  and  the  time  length 
of  the  input  signal  should  be  much  greater  than  the  period  of  the  system. 

The  procedure  of  estimating  RTFs  of  LTP  systems  using  the  first  method  is  described  in  detail.  We 
need  three  sets  of  data  to  estimate  the  transfer  functions  for  LTP  systems,  specifically  the  input  u,  output 
y,  and  time  measurements  (at  which  u  and  y  occur)  ip.  The  data  is  assumed  to  have  been  digitally 
acquired,  hence  all  the  data  points  can  be  assembled  in  a  vector  of  length  n,  where  n  is  the  total  number  of 
data  points.  Each  data  point  is  essentially  a  digital  record  of  the  measurement  of  an  analog  signal 
determined  at  some  fixed  sampling  frequency.  The  input  data  is  therefore  expressed  as 

m  =  [mj  U2  M3  •••  M„]  (14) 


y  and  i/i  are  similarly  defined. 

The  total  number  of  transfer  functions  of  the  system  that  need  to  be  identified,  ith,  have  to  be  specified 
in  advance.  An  rihxn  matrix  U  is  defined,  where  each  row  consists  of  an  appropriately  modulated  and 
Fourier  transformed  data  vector  u,  given  as 


(15) 


where  m=(  n/,-l)/2.  In  more  compact  notation,  the  U  matrix  is  really 

=  ^u^{co-mco^^  •••  u^[co)  •••  {co  + 


(16) 


where  u(m)  is  the  discrete  Fourier  transform  of  u.  Similarly,  we  define  a  Y  matrix  as  the  discrete  Fourier 
transform  of  vector  y 

Y=^{[Ti  T2  Ts  •••  T„]}  (17) 

When  identifying  LTI  systems,  we  usually  compute  the  cross-  and  power  spectral  densities  of  input  and 
output  signals,  and  estimate  the  transfer  function  of  the  system.  We  employ  similar  concepts  for  LTP 
systems,  in  that  we  define  the  power  spectral  density  of  input  signal  <I>uu,  and  the  cross  spectral  density  of 
input  and  output  signal  <I>uy  as 

Ouu  =  U*^U,Ouy=U*^Y  (18) 

where  is  complex  conjugate  transpose  of  U.  If  we  want  to  reduce  the  effect  of  noise,  we  may  use  A 
different  input  and  output  signal,  and  average  the  power-  and  cross-spectrum  as 
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(19) 


1  TV  -.TV 

k=l  k=l 


It  is  important  to  note  that  ordinary  matrix  multiplication  is  not  performed  in  the  computation  of  Ouu  and 
<I>uY  in  Equation  18.  Rather,  each  column  vector  of  U  is  transposed  to  form  a  row  vector,  and  multiplied 
with  its  corresponding  column  vector  to  yield  an  matrix.  Each  of  these  matrices  is  assembled  in  a 
three  dimensional  nhXUhxn  array  to  form  <I>uu  •  is  similarly  constructed  as  an  rihxlxn  array.  For 
example,  the  ithxrih  matrix  in  the  rth  position  of  the  third  dimension  of  Ouu  is  given  as  (denoted  as  Of/w) 


^UUi 


u-(co  -  mco^^ 
u.{a)) 
u-{w  + 


u-{w-mw^  •••  u-{w)  •••  u-ico  +  mcOp^ 


Similarly  for  the  Ouy  array,  the  rth  vector  in  its  third  dimension  is 

^^uY,=^u*{(jo-m(jOp^y.[(jo)  •••  u*[(i>)y.[(i>)  •••  u*{(jo  + m(jOp^y.[(jo)^ 

We  can  now  write  an  expression  that  represents  HTFs  of  FTP  systems  as 


^uy=^uu*G 


(20) 


(21) 


(22) 


Here,  Gis  an  rihxlxn  array  representing  HTFs  of  LTP  systems.  The  /th  vector  of  Gin  its  third  dimension 
is  given  as 

•••  Si,  So,  S-i,  •••  (23) 

Again,  the  ordinary  matrix  multiplication  is  not  performed  in  Equation  22.  The  meaning  of  the 
computation  in  Equation  22  is  that  we  can  apply  the  ordinary  matrix  multiplication  to  the  zth  matrix  (or, 

vector)  of  Ouy,  and  G .  That  is, 

^UY,  ~  ^£/£/,  g;  (24) 


It  should  be  noted  that  Equation  24  may  be  underdetermined  or  determinate,  depending  on  whether  «/,=! 
or  not.  Consider  first  the  case  «/,=!.  In  this  case,  only  the  fundamental  transfer  function  (go)  needs  to  be 

identified,  and  Ouu  and  Ouy  reduce  to  1  xl x«  arrays,  where  is  ,  and  <Puri  is  u*  {a))y,  (m). 

Each  element  of  go  is  given  as 


u*{co)y.{co) 

u,{a)f 


(25) 


which  is  the  same  relation  as  that  used  to  estimate  the  transfer  function  for  LTI  systems.  Thus,  this 
scheme  reduces  to  LTI  system  identification  if  only  one  transfer  function  is  considered  in  the  calculations. 

For  rih  >1,  however,  ^uui  matrix  is  not  positive  definite,  and  therefore  not  invertible.  That  is,  the 
problem  becomes  underdetermined,  which  implies  that  we  cannot  compute  g,  just  by  inverting  ^uui  and 
multiplying  it  to  <^uYi  in  Equation  24.  As  mentioned  earlier,  the  reason  for  the  non-invertibility  of  ^uui  is 
that  we  have  unknowns,  but  only  one  equation. 
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The  solution  we  adopt  to  overcome  this  difficulty  is  to  assume  that  the  transfer  functions  are  smooth, 
and  to  formulate  the  original  problem  as  the  minimization  of  a  cost  function  that  penalizes  the  quadratic 
estimation  error  and  the  curvature  of  RTFs.  To  describe  our  algorithm  in  more  detail  and  clearly,  the 
issues  of  the  estimation  error  and  the  curvature  of  the  transfer  functions  are  presented.  The  estimation 
error  occurs  because  we  have  made  a  simplification  in  our  analysis,  and  have  considered  only  a  few 
harmonics,  although  there  are  in  general  infinite  harmonics  in  LTP  systems.  For  instance,  suppose  a  given 
system  has  Nh  transfer  functions  of  relatively  significant  magnitudes,  but  only  nu  of  them  are  evaluated 
through  this  method.  In  that  case,  the  system  has  been  modeled  with  nu  number  of  transfer  functions,  and 
its  output  response  due  to  an  input  may  be  expressed  as 

m  _ 

Y=  (26) 

k=-m  m<|/|sM 

modeled  part  unmodeled  part 

where  m  =  (n/,-l)/2  and  M=  {Nh-l)/2.  The  unmodeled  part  essentially  appears  as  an  error,  e,  in  our  output 
equation  for  the  identification  scheme,  since  only  the  modeled  part  is  considered  in  the  identification 
procedure.  Therefore, 

m 

Y=  +  e=U^G  +  e  (27) 

k=-m 


In  addition  to  reducing  the  estimation  error  e,  we  also  need  to  make  assumptions  so  that  the 
identification  problem  is  well-posed.  In  this  regard,  we  assume  that  the  transfer  functions  are  smooth,  i.e., 
there  are  no  rapid  variations  (with  frequency)  in  the  transfer  functions.  This  is  a  reasonable  assumption, 
since  rapid  variations  with  frequency  in  a  transfer  function  usually  are  not  physical. 

In  order  to  reduce  e,  and  to  apply  the  assumption  we  have  just  made,  we  formulate  our  problem  as  the 
minimization  of  a  cost  function,  J,  that  penalizes  the  quadratic  error  and  the  curvature  of  the  transfer 
functions,  so  that  [3] 


7= min 


(Y 


(28) 


where  is  the  second  difference  operator  approximating  the  second  derivative,  given  as 


D"  = 


-2 

1 


1 

-2 


(29) 


and  a  is  a  constant.  It  is  evident  from  Equation  28  that  by  selecting  various  values  for  a,  we  can  weight 
the  second  derivatives  of  G  more  or  less  in  the  cost  function,  and  thus  penalize  the  curvature  so  that  rapid 
variations  with  frequency  in  G  are  reduced.  In  order  to  get  a  closed  form  expression  for  G,  Equation  28 
is  evaluated  by  taking  the  derivative  of  J  with  respect  to  G,  setting  it  equal  to  zero,  and  then  solving  for 
G ,  which  gives 

G=[U^U+aD"]' U^Y  (30) 

where  D  =  D  ‘D  .  Recall  that  we  had  assembled  U  U  and  U  Y  in  three  dimensional  arrays,  <I>uu  and  Ouy, 
respectively.  However,  in  order  to  numerically  evaluate  Equation  30,  we  need  to  reduce  these  arrays  into 
two-dimensional  matrices.  The  <I>uu  array  is  therefore  transformed  into  a  square,  block  diagonal,  n^xn^ 
matrix  Ouus,  in  which  the  zth  block  is  <^uui,  and  n®  =  tihxn.  Thus,  Ouus  has  the  form  as  illustrated  in  Figure 
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13(a).  Similarly,  Ouy  is  transformed  into  a  column  vector  $uy,  of  length  with  the  form  as  shown  in 
Figure  13(b). 


h«  „  H\ 


Ot/T/: 


(a)  (b) 

Figure  13:  (a)  Structure  of  Ouus  •  (b)  Structure  of  Ouy. 

The  G  matrix  is  also  transformed  into  a  vector  g,  where  g  has  n  sub-vectors  of  length  «/,,  each  of  which 
is  a  column  vector  of  G.  Equation  30  then  becomes 

g  =  [^UUs  +  «D^]"'‘I’uy  (31) 

where  indicates  the  change  of  the  dimensions  of  similar  to  the  change  in  dimensions  going  from 

U^U  to  <I>uus-  Now,  the  original  underdetermined  problem  in  Equation  24  is  transformed  into  the 
determinate  one  in  Equation  3 1  by  adopting  the  assumption  about  the  smoothness  of  the  transfer 
functions.  Using  Equation  31,  we  can  estimate  the  RTFs  of  LTP  systems,  given  the  input  signal  and 
measured  output  response.  Here,  the  weighting  factor  a  should  be  selected  carefully  such  that  the  transfer 
functions  become  smoothed  out  enough,  while  capturing  the  important  features  of  the  systems.  Note  that 
too  much  smoothing,  i.e.,  too  high  value  of  a,  may  overestimate  the  system  damping,  and  therefore,  the 
results  of  the  system  identification  may  be  inaccurate. 

The  advantage  of  the  first  method  is  that  we  can  generate  an  input  signal  quite  easily  to  estimate  the 
RTFs.  That  is,  the  results  of  system  identification  using  this  method  are  quite  insensitive  to  the  input 
signal,  as  long  as  the  input  signal  has  enough  frequency  components  and  time  length.  As  will  be  shown  in 
the  next  section,  this  is  a  significant  advantage  over  the  second  method,  whose  solutions  are  sensitive  to 
small  errors  in  the  input  signal. 

However,  the  first  approach  needs  some  ad-hoc  manners  to  transform  the  originally  underdetermined 
problem  into  the  determinate  one,  and  therefore  to  estimate  the  RTFs.  In  addition,  this  approach  requires 
much  computer  memory  and  intensive  computational  procedure,  because  we  should  take  into  account  the 
power  spectral  density  <I»uus  and  cross-spectral  density  matrix  Ouy  of  input  and  output  spectrum  with  all 
frequencies  simultaneously. 


Second  method 

The  second  approach  we  developed  to  solve  the  underdetermined  problem  is  to  apply  inputs  and 
measure  the  resulting  outputs,  in  order  to  estimate  «/,  harmonic  transfer  functions.  We  can  thus  form  hh 
independent  equations  of  the  form  of  Equation  12,  and  will  subsequently  be  able  to  evaluate  the  desired 
transfer  functions.  This  approach  allows  us  to  avoid  the  use  of  the  heuristic  assumption  that  the  transfer 
functions  are  smooth,  and  therefore,  we  can  address  the  problem  more  naturally.  In  order  to  use  this 
approach  successfully,  however,  we  should  be  extremely  careful  in  designing  input  signals  to  excite  LTP 
systems,  compared  with  the  first  approach.  That  is,  it  is  very  important  to  take  into  consideration  the  time 
of  application  of  each  input  relative  to  the  system  period.  This  is  due  to  the  time-varying  nature  of  the 
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system  dynamics  during  one  period.  If  the  system  behavior  is  to  be  completely  characterized,  the  system 
needs  to  be  excited  with  appropriate  input  signals  at  various  times  in  its  period  T.  When  «/,  transfer 
functions  are  being  evaluated,  at  least  the  same  number  of  identical  input  signals  should  be  applied  that 
are  evenly  spread  out  over  the  system  period.  Figure  14  depicts  when  each  of  the  «/,  inputs  should  be 
initiated,  relative  to  the  system  period. 


'  Beginning  of 
'  ^  a  new  period 
of  the  system 


Beginning  of 
a  new  period 
of  the  system 


Figure  14:  Input  signals  initiated  at  appropriate  time  intervals  over  the  system  periods. 


We  define  7^,  as  the  amount  of  time  (in  seconds)  elapsed  between  the  beginning  of  a  new  period  of 
the  system  and  the  start  of  the  first  input  signal  since  the  beginning  of  that  particular  period.  In  this 
specific  case,  where  we  have  just  established  that  the  «/,  inputs  should  be  spaced  uniformly  apart,  Td  is 
given  as 


III 


(32) 


The  first  input  should  have  zero  delay  between  the  start  of  a  system  period  and  its  time  of  initiation.  The 
response  of  the  LTP  system  output,  >>0(05  due  to  the  first  input,  u{t),  is  given  by 


m 

>^0(0= 


k=-m 


(33) 


where  m  =  («/,-l)/2,  denotes  the  convolution  operator,  and  gk(t)  is  the  kth  harmonic  impulse  response  of 
the  system.  The  input  U=F{u(t)},  and  output  Ya=F{yQ{t)},  to  the  system  for  n/,=3  in  the  frequency  domain 
are  modeled  as  depicted  in  Figure  12. 

For  the  second  chirp,  there  should  be  a  delay  of  Td  seconds  between  the  start  of  a  system  period  and 
its  time  of  application,  so  the  actual  input  will  be  u(t-Td).  The  system  output  due  to  this  second  chirp,  yi(t), 
is  given  by 

m 

yi{^-Td)=  (34) 


Repeating  the  procedure  described  above,  the  system  output  due  to  the  «th  chirp,  y„(0,  will  be  expressed 
as  in  the  time  domain 


T, 


m 

k=-m 


<n<n^ 


(35) 


Therefore,  the  relation  between  the  nth  chirp  input  and  its  resulting  system  output  in  the  frequency 
domain  is  given  as 


m 

=  Jg, (>)(/(>  -  0  <n< 


n,. 


(36) 


k=-m 


24 


Table  of  Contents  (including  Appendices) 


TABLE  OF  CONTENTS  (INCLUDING  APPENDICES) . 2 

TABLE  OF  ILLUSTRATIONS  AND  TABLES . 3 

1.  STATEMENT  OF  PROBLEM  STUDIED . 4 

2.  SUMMARY  OF  IMPORTANT  RESULTS . 4 

3.  PAPERS  SUBMITTED  OR  PUBLISHED . 7 

4.  SCIENTIFIC  PERSONNEL . 7 

5.  REPORT  OF  INVENTIONS . 7 

6.  BIBLIOGRAPHY  (CITED  ABOVE  AND  IN  APPENDICES) . 8 

A  FLAP  POSITION  CONTROL  LOOP . 9 

A.  1  Basic  Concept  of  Flap  Position  Control  Loop . 9 

A.2  Transfer  Function  of  Actuator  Systems . 9 

A.3  Design  of  Controllers  and  Simulation  Results . 10 

A. 4  Discussions . 16 

B  SYSTEM  IDENTIFICATION  OF  LINEAR  TIME  PERIODIC  SYSTEMS . 16 

B .  1  Frequency  Response  of  LTP  Systems . 17 

B.2  System  Identification  Schemes  for  LTP  systems . 19 

B .  3  V ALID ATION  OF  SYSTEM  IDENTIFICATION  METHODS . 27 

C  SUMMARY  OF  WIND  TUNNEL  TEST  AT  LANGLEY . 34 

C .  1  System  Identification  of  Rotor  Systems . 35 

C.2  Closed-loop  Experimental  Results . 39 


2 


where  Y„(Jw),  U(Jw),  and  Gk(J(o)  are  Fourier  transform  of y„{t),  u{t),  &nd,gk{f),  respectively.  It  should  be 
noted  that  the  Fourier  transform  of  u{t  -  nT^  )e^^”'’^is  jo)  -  .  Combining  all  of  the 

system  response,  the  output  vector  Y^=[Fo  Yi  ... Y„h.i]  can  be  expressed  as 


•  u,  •• 

..  . 
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•  U,  •• 
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u  w 


n(nj-l) 


Gk 
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(37) 


jcO  T  ^  'T' 

where  U„=U{ja)  +jm(Op),  and  W  =  e  "  .  Multiplying  U  to  both  sides  of  Equation  37  yields 

(38) 

where  Ouy  =  U  ^  Y  and  Ouu  =  U  ^  U.  As  in  the  first  method,  we  can  use  N  different  input  and  output 
signal  and  average  the  power-  and  cross-spectrum,  in  order  to  reduce  the  effect  of  noise.  In  that  case,  <I>uu 
and  OuY  are  given  as 


k=l  k=l 


(39) 


In  contrast  to  the  first  approach,  the  power  spectral  density  of  the  input  signal  <I>uu  is  not  singular  in  the 
second  approach.  The  reason  for  the  invertibility  of  <I>uu  is  that  we  consider  at  least  the  same  number  of 
the  input  signals  as  that  of  the  HTFs  to  be  identified.  Therefore,  we  can  compute  the  HTFs  vector  G  from 
Equation  38  just  by  inverting  Ouu  and  multiplying  it  to  Ouy,  as 

G  =  OuuOuy  (40) 


Equation  40  can  be  further  simplified,  so  that  the  kYa  HTF,  Gk,  can  be  obtained  after  some  mathematical 
manipulations,  as  (for  A=l) 


U*[jm  -  jkcOp )  2  Y„  {ja))W 


kn 


«=0 


u[i(o-ika>p) 


(41) 


It  should  be  noted  that  the  computation  of  G  in  Equation  40  can  be  performed  at  each  frequency 
independently.  That  is,  Ouu  is  an  nh'x-nh  matrix,  and  therefore,  the  number  of  data  points  is  not  involved  in 
the  inverse  matrix  computation  of  <5uu-  This  makes  a  critical  difference  between  the  first  and  the  second 
method.  Recalling  from  Equation  31  is  an  nnhxnuh  matrix  (n  is  the  number  of  data  points),  we  have 
observed  that  the  first  approach  requires  us  to  consider  the  whole  frequency  response  simultaneously  to 
estimate  HTFs. 

Equation  40  will  yield  HTFs  of  LTP  systems  with  reasonable  accuracy,  if  the  effect  of  noise  is  not 
significant.  However,  if  the  system  output  is  so  severely  corrupted  by  noise  that  the  averaging  in  Equation 
39  is  not  enough  to  reduce  the  effect  of  noise,  some  smoothing  schemes  should  be  employed  to  generate 
transfer  functions  with  acceptable  accuracy.  First,  we  can  convolve  the  averaged  power-  and  cross¬ 
spectrum  in  the  frequency  domain  with  smoothing  window  functions.  The  smoothed  power-  and  cross¬ 
spectrum  are  given  by 

=  ^vv  *  (42) 
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where  w((o)  is  a  window  function  and  is  a  convolution  operator,  both  defined  in  the  frequency  domain. 
Various  window  functions  can  be  used  here,  such  as  a  rectangular,  Hanning,  Balett  window,  etc., 
depending  on  the  applications.  The  smoothed  HTF  is  obtained  by  taking  the  ratio  of  the  smoothed  cross- 
and  power  spectrum  as 

G  =  Ouu4’uy  (43) 

The  smoothing  technique  in  Equation  43  reduces  the  effect  of  measurement  noise  quite  effectively  and 
leads  to  an  accurate  estimate  of  the  transfer  functions,  as  long  as  the  smoothing  window  function  is 
selected  carefully.  On  the  other  hand,  the  smoothing  technique  in  Equation  43  may  not  be  enough  if  the 
input  signal  is  also  corrupted  by  noise,  or  has  some  problems  in  terms  of  frequency  components  and 
synchronization.  The  corruption  of  the  input  signal  may  occur  when  there  is  any  unexpected  problem  in 
the  electronics  hardware  or  frequency  analyzer  to  generate  the  input  signal.  Also,  the  results  of  the  system 
identification  may  be  severely  degraded,  if  the  input  signal  does  not  have  enough  frequency  components 
and/or  the  application  of  each  input  signal  is  not  exactly  synchronized  with  the  system  period  or  sampling 
time  of  the  DSP  system.  As  evident  from  Equation  41,  we  should  consider  the  frequency  range  of 
u{ja)  -  jkojp  ),  the  shifted  Fourier  transform  of  the  input  signal,  as  well  as  U(ju)),  in  order  to  estimate  the 

HTFs  with  good  precision.  If  the  input  signal  does  not  have  enough  frequency  components, 

{/(/to  -  jkcjp  )  may  have  very  low  magnitude  in  the  bandwidth  of  interest,  while  U(J(n)  has  reasonable 

magnitude,  so  that  the  estimate  of  the  kXh  HTF  Gt  may  be  very  inaccurate. 

In  order  to  reduce  the  effect  of  the  inaccuracy  in  the  input  signal  on  the  system  identification 
performance,  we  devised  another  smoothing  technique,  in  which  some  small  penalty  is  added  to  the 
power  spectral  density  of  the  input  signal  before  inverting  the  matrix.  In  other  words, 

6-('i>™  +  £i)"''S'uv  (44) 

Here,  I  is  an  identity  matrix  compatible  with  and  £  is  a  constant.  This  smoothing  process  reduces 
the  magnitudes  of  the  HTF  a  little  bit,  but  reduces  noise  in  the  HTFs  quite  significantly.  The  role  of  the 
small  penalty  e  is  to  make  the  estimate  of  the  transfer  functions  robust  to  the  error  or  inaccuracy  in  the 
input  signal. 

As  will  be  shown  in  the  next  section,  noise  in  the  input/output  signal  corrupts  the  estimate  of  the 
transfer  functions  quite  significantly,  and  makes  it  difficult  to  distinguish  between  the  HTFs  with 
negligible  magnitude  but  much  noise,  and  those  with  much  magnitude  but  significant  variations. 

However,  we  can  distinguish  them  easily  using  these  smoothing  processes  in  Equation  43  and  44.  That  is, 
HTFs  insensitive  to  the  window  function  w(co)  and/or  small  penalty  £  are  real,  while  those  sensitive  to  the 
smoothing  process  are  artifact. 

When  employing  the  second  approach  to  estimate  HTFs,  it  should  be  pointed  out  that  the  output 
measurement  due  to  each  input  must  be  obtained  by  allowing  the  response  to  settle  down  significantly 
before  the  next  input  signal  is  initiated.  In  that  case,  it  can  be  reasonably  assumed  that  To  is  due  to  the  first 
input,  Ti  due  to  the  second  input,  and  so  on. 

In  contrast  to  the  first  method,  the  second  method  doesn't  require  much  computer  memory  or 
intensive  computation  work,  because  we  can  compute  the  transfer  functions  at  each  frequency 
independently,  without  considering  the  whole  frequency  response  simultaneously.  On  the  other  hand,  in 
the  first  method,  we  should  generate  several  'big'  matrices,  which  cause  much  memory  and  long 
computation  time,  because  the  ad-hoc  manner  in  the  first  method  requires  that  we  should  take  into 
account  the  power  spectral  density  and  cross-spectral  density  matrix  of  input  and  output  spectrum  with  all 
frequencies  simultaneously. 
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B.3  Validation  of  System  Identification  Methods 

Once  we  have  developed  the  methodologies  to  identify  LTP  systems,  it  is  important  to  validate  our 
identification  schemes  before  applying  them  to  actual  systems.  For  this  reason,  we  selected  a  simple  LTP 
system  model,  a  Mathieu  equation,  whose  harmonic  transfer  functions  can  be  found  analytically,  and 
verified  the  accuracy  of  our  schemes  using  the  model.  We  first  analyze  the  model  by  computing  their 
HTFs  analytically.  Then,  we  apply  our  identification  schemes  to  estimate  the  HTFs  of  the  model  using  the 
designed  input  data  and  simulated  output  data,  and  assuming  that  the  model  is  given  as  a  'blackbox'. 
Finally,  we  compare  the  analytic  solutions  with  the  results  obtained  using  our  identification  schemes. 


Theoretical  Analysis 


The  Mathieu  Equation  is  a  well-known  equation  representing  an  LTP  system.  The  canonical  form  of 
the  lossy  Mathieu  equation  that  we  employed  for  our  analysis  is 

x{t)  +  2^co^x(t)  +  {l  +  2l3cosco^t^colx(t)=0  (45) 


Note  that  the  parameter  (3  in  the  equation  determines  the  importance  of  the  periodic  effects.  That  is,  the 
system  is  essentially  LTI  for  small  f},  while  the  periodic  effects  are  significant  for  large  fi. 

In  state-space,  the  system  matrices  are 


A(0  = 


0  1 
-{l  + 2/3  cos  (Optical  -2^co^ 

The  harmonics  of  the  dynamics  matrix,  A(t),  are 


,B(0  = 


,C(0  =  [1  1] 
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Since  there  are  no  periodic  elements  in  B(t)  and  C(t)  matrices.  Bo  and  Co  are  the  same  as  the  matrices 
B(t)  and  C(t)  given  in  Equations  46,  while  all  the  higher  harmonics  are  zero.  Using  Equation  1 1,  we  can 
calculate  the  HTFs  of  the  Mathieu  equation.  Figure  15  shows  four  HTFs  of  the  Mathieu  equation.  Other 
HTFs  have  negligible  magnitude,  so  they  are  not  shown  in  Figure  15.  The  system  parameters  used  in  the 
simulation  are  (3=0.3,  ^=0.01,  a)„=20n  rad/sec  (10  Hz),  and  a)p=S0n  rad/sec  (40  Hz). 

The  HTFs  shown  in  Figure  15  show  the  importance  of  LTP  system  analysis  and  identification.  If  the 
required  bandwidth  for  this  system  is  below  20  Hz,  we  can  assume  this  system  is  LTI  with  the  transfer 
function  Go,  and  apply  various  LTI  control  design  methods  to  achieve  the  performance  requirement, 
because  Go  has  much  higher  magnitude  than  other  transfer  functions  at  least  by  20  dB  below  20  Hz. 
However,  if  the  bandwidth  of  interest  is  higher  than  20  Hz,  the  LTI  assumption  does  not  hold,  and  we 
should  analyze  and  design  this  system  as  LTP  systems,  because  other  transfer  functions,  such  as  Gi  and 
G.i,  have  similar  or  even  higher  magnitude  than  Go. 
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Figure  15:  Frequency  responses  of  the  lossy  Mathieu  equation. 

Empirical  Analysis 

After  we  analyzed  the  Mathieu  equation  model  analytically,  we  used  our  identification  algorithms  to 
estimate  the  HTFs  of  the  system.  In  order  to  evaluate  the  identification  methods  thoroughly  and  make  the 
system  identification  procedure  look  more  realistic,  we  considered  three  different  cases.  In  the  first  case, 
we  assume  that  both  the  input  signal  and  output  measurement  are  perfect,  and  there  is  no  noise  at  all.  By 
doing  that,  we  can  make  sure  that  the  LTP  system  can  be  actually  identified  using  our  algorithms,  and  we 
can  analyze  the  effects  of  design  parameters  on  the  accuracy  of  the  algorithms.  On  the  other  hand,  we 
assume  that  the  output  measurement  is  corrupted  by  noise  in  the  second  case.  We  artificially  generate  the 
white  noise,  and  add  it  to  the  output  signal  of  the  Mathieu  equation.  In  the  third  case,  the  input  signal  as 
well  as  the  output  signal  are  assumed  to  be  corrupted  by  noise.  We  used  the  actual  input  signal  employed 
used  to  estimate  the  HTFs  of  a  helicopter  rotor  during  the  wind  tunnel  test  at  Langley.  By  using  these 
artificial  noise  in  input  and/or  output  signal,  combined  with  the  model,  we  can  see  whether  our  algorithms 
can  be  used  to  identify  the  actual  physical  systems.  The  following  sections  present  the  results  of  the 
analysis  for  these  three  cases. 

(1)  Perfect  Measurement  and  Perfect  Input  Signal 

First,  it  was  assumed  that  both  the  input  and  output  signal  are  perfect,  and  the  two  identification 
algorithms  were  used  to  estimate  the  HTFs  of  the  Mathieu  equation.  In  order  to  apply  the  algorithms 
successfully,  especially  the  second  method,  the  input  signal  should  be  designed  carefully.  Considering 
that  we  want  to  estimate  four  HTFs  in  this  example,  we  need  to  use  at  least  four  chirp  signals  to  estimate 
HTFs  with  good  precision.  In  this  example,  we  chose  five  chirp  signals,  each  of  which  lasts  10  seconds 
and  sweeps  sine  waves  from  0  to  200  Hz.  Each  chirp  in  the  input  sequence  differed  in  phase,  by  equal 
amounts  over  a  period.  Since  we  use  five  chirps,  the  phase  difference  between  two  consecutive  signals  is 
271/5.  We  also  add  enough  zero  input  to  each  chirp  signal,  while  taking  into  account  the  relative  phase,  to 
allow  the  response  to  decay  enough  before  the  next  input  signal  is  initiated. 

Figure  16  shows  the  input  signals  designed  for  this  example,  and  the  resulting  output  response  of  the 
system.  The  same  parameters  as  in  Section  B.3  were  used  to  simulate  the  Mathieu  equation  model.  As 
expected,  the  output  response  corresponding  to  each  of  the  chirp  signals,  which  differ  only  in  phase,  is 
slightly  different  in  magnitude. 
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Figure  16:  Simulated  input  signals  and  corresponding  response  of  the  Mathieu  equation. 

We  now  apply  our  identification  schemes  to  the  input  and  output  response  in  Figure  16,  to  estimate 
four  HTFs.  Figure  17(a)  and  (b)  show  the  HTFs  of  the  Mathieu  equation  model  estimated  using  the  first 
and  the  second  method,  respectively,  along  with  the  analytic  transfer  functions.  When  applying  the  first 
method,  the  input  signal  in  Figure  16  was  treated  as  one  long  input  signal,  and  the  least  square  technique 
was  applied  to  yield  a  well-defined  identification  problem.  Also,  the  fourth  difference  operator  was  used 
to  smooth  the  responses,  and  various  values  of  the  weighting  factor  a  were  tested  to  investigate  the  effect 
of  a  on  the  identification  results  and  find  a  value  resulting  in  the  acceptable  accuracy  in  the  transfer 
functions.  The  result  in  Figure  17(a)  corresponds  to  a  =100,  but  we  found  that  the  results  of  the  first 
identification  scheme  are  quite  insensitive  to  a.  That  is,  the  transfer  functions  remain  almost  unchanged 
when  a  varies  from  10“’  to  10*. 

On  the  other  hand,  the  second  method  treats  the  five  chirps  in  Figure  16  as  five  different  inputs  with 
different  phase,  and  computes  HTFs  by  inverting  directly  the  power  spectral  density  of  the  input 
spectrum.  The  HTFs  in  Figure  17(b)  were  obtained  using  Equation  43,  i.e.,  without  using  any  smoothing 
scheme,  because  it  was  assumed  that  there  is  no  noise  in  the  input  and  measurement  signal.  Also,  we 
could  use  only  five  chirps  in  the  input  signal  and  estimate  four  HTFs  quite  successfully,  but  we  should 
consider  in  general  as  many  chirp  signals  as  possible  to  improve  the  accuracy  of  the  identification 
algorithms. 

Figure  17(a)-(b)  show  that  the  analytic  transfer  function  matches  well  with  estimated  transfer 
functions.  The  discrepancies  between  analytic  and  estimated  transfer  functions  are  so  small  that  their 
effect  will  be  negligible. 
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Figure  17:  Estimated  harmonic  transfer  functions  of  the  Mathien  eqnation  system  obtained  nsing  (a)  the  first 
method,  (b)  the  second  method.  Solid  and  dotted  line  corresponds  to  analytical  and  empirical  response, 
respectively. 


(2)  Imperfect  Measurement  and  Perfect  Input  Signal 

In  the  second  case,  we  assumed  that  the  output  measurement  is  corrupted  by  noise,  while  the  input 
signal  is  carefully  designed  and  noiseless.  In  order  to  mimic  the  noisy  measurement,  we  artificially 
generated  the  white  noise,  and  added  it  to  the  output  signal  of  the  Mathieu  equation.  The  main  motivation 
of  this  section  is  to  investigate  the  effects  of  the  sensor  noise  on  the  performance  of  our  identification 
schemes,  and  figure  out  how  the  design  parameters  can  reduce  the  noise  effect. 

Figure  18  shows  the  input  signal  designed  for  this  example,  which  is  the  same  as  that  in  the  previous 
section,  and  the  resulting  output  response  with  the  noise.  By  comparing  with  the  results  in  Figure  16,  we 
can  see  clearly  that  the  output  response  in  Figure  18  is  severely  corrupted  with  the  noise.  The  two  FTP 
system  identification  methods  were  then  applied  to  the  input  signal  and  output  response  in  Figure  18.  Ten 
input  and  output  signals  were  averaged  in  the  frequency  domain  using  the  cross-  and  power  spectrum,  to 
reduce  the  effects  of  the  measurement  noise.  The  results  of  applying  the  first  identification  method  are 
shown  in  Figure  19,  with  various  weight  factors.  The  figure  shows  clearly  the  effect  of  the  weighting 
factor  on  the  estimate  of  RTFs,  when  they  are  corrupted  with  noise.  It  is  evident  from  the  figure  that  a 
high  weighting  factor  smoothes  out  the  RTFs,  by  the  sacrificing  rapid  variations  with  frequency  or 
overestimating  the  system  damping.  When  a<10*  (Figure  19(a)-(c)),  there  is  little  change  in  Go(s)  and  Gi 
(s),  which  have  much  higher  magnitudes  than  other  transfer  functions.  Rowever,  when  a  =10^  (Figure 
1 9(d)),  we  can  see  that  the  transfer  functions  are  too  smoothed  out,  that  is,  the  system  damping  is 
overestimated.  This  fact  implies  that  we  should  use  engineering  judgment  in  selecting  the  weight  factor, 
so  that  the  transfer  functions  are  smoothed  enough,  while  preserving  their  important  features. 
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Figure  18:  Simulated  input  signals  and  corresponding  response  of  the  Mathieu  equation.  The  output  is 
corrupted  by  white  noise. 


Figure  19:  Harmonic  transfer  functions  of  the  Mathieu  equation  system  obtained  using  the  first  method  with 
various  a.  (a)  a  =1,  (b)  a  =10^  (c)  a  =10*,  (d)  a  =10®. 
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Figure  20:  Harmonic  transfer  functions  of  the  Mathieu  equation  system  obtained  using  the  second  method. 

(a)  Equation  40  (without  smoothing  technique),  (b)  Equation  43  (with  the  first  smoothing  technique),  (c) 
Equation  44  (with  the  first  and  second  smoothing  techniques). 

Figure  20  shows  the  results  of  using  the  second  identification  method  to  estimate  the  HTFs,  without 
any  smoothing  (Figure  20(a)),  with  smoothing  in  Equation  43  (Figure  20(b)),  and  with  smoothing  in 
Equation  44  (Figure  20(c)).  We  can  see  that  the  smoothing  technique  in  Equation  43  significantly  reduces 
the  effect  of  measurement  noise,  and  improves  the  accuracy  of  the  identification  method.  On  the  other 
hand,  it  is  difficult  to  see  the  difference  between  the  results  in  Figure  20(b)  and  Figure  20(c),  that  is, 
smoothing  technique  in  Equation  43  and  44,  because  it  is  assumed  that  the  input  signal  is  well  designed 
without  any  noise  in  this  section.  As  mentioned  earlier,  the  smoothing  scheme  in  Equation  44  was  devised 
to  reduce  the  effect  of  any  inaccuracy  in  the  input  signal  on  the  performance  of  the  identification  method. 
The  advantage  of  using  the  smoothing  in  Equation  44  will  be  evident  in  the  next  section. 

(3)  Imperfect  Measurement  and  Imperfect  Input  Signal 

In  the  previous  section,  it  was  assumed  that  the  measurement  is  noisy  while  the  input  signal  is  perfect, 
and  the  effects  of  the  measurement  noise  on  the  performance  of  our  identification  schemes  were 
investigated.  In  this  section,  we  assume  that  the  measurement  is  corrupted  with  noisy,  and,  in  addition,  the 
input  signal  is  faulty,  in  that  it  is  also  noisy  and  does  not  have  enough  frequency  components.  We  used 
the  actual  input  signal  used  during  the  wind  tunnel  test  at  Langley  to  estimate  the  harmonic  transfer 
functions  of  a  helicopter  rotor,  which  is  linearly  varying  sine-sweep  signal  from  5  Hz  to  70  Hz.  The  input 
signal,  as  well  as  the  resulting  output  response  of  the  Mathieu  equation,  are  shown  in  Figure  21.  Although 
not  clear  from  the  figure,  the  input  signal  is  corrupted  with  several  spikes,  due  to  unexpected  behavior  of 
the  experimental  facility.  As  for  the  measurement  noise,  the  same  white  noise  used  in  the  previous  section 
was  added  to  the  output  response  of  the  Mathieu  equation. 

Figure  22  shows  the  results  of  applying  the  first  identification  method  to  the  input  and  output 
response  in  Figure  21,  with  the  same  weights  as  in  the  previous  section.  Compared  with  the  results  in 
Figure  19,  we  can  see  that  the  effect  of  noise  is  more  severe  here  but  not  much,  and  the  smoothing  with 
adequate  weighting  factor  can  effectively  reduce  it.  On  the  other  hand,  the  results  in  Figure  23,  where  the 
second  identification  method  was  applied,  show  that  the  effect  of  noise  is  much  more  significant  if  the 
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input  signal  is  faulty,  compared  with  the  results  in  Figure  20.  That  is,  in  contrast  to  the  first  method,  the 
second  identification  method  is  sensitive  to  the  input  signal.  Also,  the  difference  between  the  results  in 
Figure  23(b)  (smoothing  technique  in  Equation  43)  and  Figure  23(c)  (smoothing  technique  in  Equation 
44)  is  now  apparent,  which  was  not  in  the  previous  section.  By  using  the  smoothing  in  Equation  44,  as 
well  as  in  Equation  43,  we  can  reduce  considerably  the  effect  of  inaccuracy  in  the  input  signal  on  the 
performance  of  the  second  method.  Therefore,  we  should  be  very  careful  in  designing  the  input  signal  for 
the  second  method  and  use  the  smoothing  methods  in  Equation  43  and,  especially  Equation  44,  to  make 
the  method  more  robust  to  unexpected  problems  in  the  input  signal. 


Figure  21:  Input  signal  intentionally  designed  faulty,  and  corresponding  response  of  the  Mathieu  equation. 
Both  input  and  output  signal  arc  corrupted  with  noise. 


Figure  22:  Flarmonic  transfer  functions  of  the  Mathieu  equation  system  obtained  using  the  first  method  with 
various  a.  (a)  a=l,  (b)  a=10^  (c)  a  =10*,  (d)  a  =10*. 
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Figure  23:  Harmonic  transfer  functions  of  the  Mathieu  equation  system  obtained  using  the  second  method, 
(a)  Equation  40  (without  smoothing  technique),  (b)  Equation  43  (with  the  first  smoothing  techniqne).  (c) 
Eqnation  44  (with  the  first  and  second  smoothing  techniques). 


C  Summary  of  Wind  Tunnel  Test  at  Langley 

While  Boeing  prepared  for  the  flight  test,  we  collaborated  with  Prof  Carlos  Cesnik’s  research  group 
at  University  of  Michigan  to  support  his  NASA-supported  program  for  the  Active  Twist  Rotor  (ATR)  in 
forward  flight  in  the  Transonic  Dynamics  Tunnel  (TDT)  at  NASA  Langley.  Figure  24  shows  the  ATR 
system  mounted  on  the  Aeroelastic  Rotor  Experimental  System  (ARES)  helicopter  testbed  with  in  the 
TDT. 


Eigure  24:  ATR  system  mounted  on  the  ARES  helicopter  testbed  with  in  the  TDT. 

The  collaboration  was  carried  out  to  test  our  system  identification  methodology  and  continuous-time 
higher  harmonic  control  (HHC)  algorithms  on  actual  rotor.  For  the  wind  tunnel  test,  we  implemented  the 
system  identification  algorithms  in  Simulink  to  estimate  the  transfer  function  of  the  rotor,  and  designed 
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continuous-time  HHC  algorithms  incorporating  an  anti-windup  scheme.  We  also  transferred  our  system 
identification  and  HHC  design  technology  to  Prof  Cesnik’s  group  to  support  the  wind  tunnel  test  for  the 
ATR  at  NASA  Langley.  This  section  summarizes  the  results  of  the  flight  test  at  Langley. 


C.l  System  Identification  of  Rotor  Systems 

We  applied  our  identification  schemes  described  in  Section  B  to  the  ATR,  which  consists  of  a  fully 
articulated  four  blades  incorporating  active  materials  as  actuators.  We  estimated  the  transfer  functions  of 
the  ATR  system  at  various  advance  ratio  (^=0. 14,  0.2,  0.267,  0.333),  by  applying  electric  field  to  the 
active  materials  made  of  Active  Fiber  Composite  (AFC)  packs  in  each  ATR  blade,  and  measuring  the 
system  output  from  a  force  sensor  mounted  inside  the  helicopter  hub.  In  order  to  estimate  the  HTFs  of  the 
ATR  system  using  our  identification  schemes,  we  carefully  designed  the  input  signal  to  excite  the  system 
and  measure  the  output  response,  considering  the  system  period,  7=2^1  (Op,  and  the  number  of  HTFs  to  be 
identified,  «/,.  We  used  a  linearly  varying  sine-sweep  signal  from  5  Hz  to  70  Hz  for  each  chirp  input,  since 
the  main  frequency  of  interest  in  our  experiment  is  4  /rev,  which  is  approximately  46  Hz.  The  whole  input 
signal  consists  of  20  chirp  signals,  each  of  which  lasts  20  seconds  and  has  equally-spaced  relative  phase. 
We  used  20  chirp  signals,  which  may  appear  more  than  adequate  in  this  case,  in  order  to  improve  the 
accuracy  of  the  identification  algorithms.  It  should  be  noted  that  considering  more  chirps  in  the  input 
signal  achieves  better  identification  results.  Therefore,  it  becomes  more  critical  to  consider  as  many  chirps 
as  possible  when  applying  the  algorithms  to  actual  systems,  because  the  measurement  output  may  be 
seriously  corrupted  by  sensor  noise  in  a  real  world.  Figure  25  shows  some  parts  of  the  input  signal  and  the 
corresponding  response  of  the  ATR  system. 

Before  applying  our  identification  schemes  to  estimate  the  HTFs  of  the  rotor  system  (collective,  sine 
cyclic,  cosine  cyclic)  at  various  advance  ratios,  we  selected  one  specific  case,  i.e.,  collective  transfer 
function  at  ^=0.333,  and  tested  our  two  system  ID  methods  in  detail.  We  followed  the  same  procedure  in 
Section  B.3,  in  order  to  see  whether  our  methods  yield  the  physically  meaningful  results.  In  the  first 
method,  we  treated  the  input  signal  consisting  of  20  chirps  as  one  long  input  signal,  and  employed  the 
least  square  technique  to  yield  a  well-posed  identification  problem.  Since  there  is  no  decisive  way  of 
determining  the  optimal  value  of  the  weighting  factor  a,  we  began  with  a  =10^  and  increased  the  value 
exponentially  to  find  a  reasonable  value  achieving  the  acceptable  accuracy.  Figure  26(a)-(d)  show  the  five 
collective  HTFs  at  ^-0.333,  which  were  obtained  using  the  first  method  with  a  =10®,  10^,  lO’^,  and  10^®, 
respectively.  The  case  for  a  =10^  was  not  shown  here,  because  the  resulting  transfer  functions  were  too 
noisy. 
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Figure  25:  Designed  input  signal  and  corresponding  response  of  the  ATR  system. 

Figure  26  shows  clearly  the  effect  of  the  weighting  factor  on  the  HTFs.  It  is  evident  from  the  figure 
that  a  high  weighting  factor  smoothes  out  the  HTFs,  while  sacrificing  the  rapid  variations  with  frequency 
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in  them.  Since  it  is  often  difficult  to  figure  out  whether  the  rapid  variation  is  due  to  the  noise  or  exists 
physically,  we  should  test  various  values  of  the  weighting  factor  and  compare  their  results.  We  can  also 
observe  from  Figure  26  that  there  exist  RTFs  quite  insensitive  to  the  weighting  factor.  That  is,  we  can  see 
that  Go(s),  Gi(s),  and  G2(s)  do  not  change  much  when  we  vary  the  weighting  factor  a  from  10®  to  lO'®, 
while  G.2(s)  and  G.i(s)  change  significantly  depending  on  a.  When  a  =10®,  the  magnitude  of  G.2(s)  and 
G.i(s)  is  similar  to,  or  even  higher  than,  that  of  the  other  three  RTFs,  but  G.2(s)  and  G.i(5)  become 
negligible  as  a  increases.  What  this  fact  may  imply  is  that  Go(s),  Gi(s),  and  G2(s)  (which  do  not  depend 
much  on  a)  are  the  transfer  functions  that  actually  exist  in  the  system,  while  G.2(s)  and  G.i(s)  (which  are 
very  sensitive  to  a)  do  not  exist  physically,  or  have  so  low  magnitudes  that  they  cannot  be  identified  or 
can  be  neglected. 

In  contrast  to  the  first  method,  the  second  method  treats  the  20  chirps  (some  of  them  are  shown  in 
Figure  25)  as  20  different  inputs  with  different  phase,  and  measures  20  different  output  response  to 
estimate  RTFs.  Figure  27(a)  shows  the  five  RTFs  obtained  using  Equation  40,  that  is,  without  using  any 
smoothing  scheme.  Figure  27(b)  shows  the  five  RTFs  obtained  using  the  smoothing  technique  in 
Equation  43,  while  Figure  27(c)  shows  the  same  results  using  Equation  44.  Figure  27  clearly  indicates 
how  the  smoothing  procedures  affect  the  resulting  RTFs  estimated  using  the  second  method.  It  can  be 
observed  from  the  figure  that  the  smoothing  procedures  in  the  second  method  reduce  the  noise  in  the 
RTFs  significantly,  and  enable  us  to  distinguish  between  the  RTFs  with  negligible  magnitude  but  much 
noise,  and  those  with  much  magnitude  but  significant  variations.  Figure  27  shows  that  the  use  of  the 
smoothing  procedures  make  little  differences  in  Go(s),  Gi(s),  and  G2(s),  while  there  is  significant  change 
in  G.2(s)  and  G.i(s)  depending  on  the  window  function  and  the  penalty  e.  This  observation  implies  that 
RTFs  insensitive  to  the  smoothing  exist  actually  in  the  system,  while  those  sensitive  to  smoothing  may  be 
an  artifact  of  inaccurate  measurement.  It  should  be  also  noted  that  we  could  draw  the  same  conclusions 
using  the  first  method. 


Figure  26:  Harmonic  transfer  functions  of  the  ATR  system  obtained  using  the  first  method  with  various  a.  (a) 
a  =10®,  (b)  a=10®,  (c)  a=10'^  (d)  a=10'^ 
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Figure  27:  Harmonic  transfer  functions  of  the  ATR  system  obtained  using  the  second  method,  (a)  Without 
smoothing  procedure,  (b)  With  smoothing  procedure. 


Figure  28:  Harmonic  transfer  functions  of  the  ATR  system,  (a)  Collective,  (b)  Cosine  cyclic,  (c)  Sine  cyclic. 
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Figure  29:  Effect  of  advance  ratio  on  the  HTFs  of  the  ATR  system,  (a)  Collective,  (b)  Cosine  cyclic,  (c)  Sine 
cyclic. 


After  we  verified  our  identification  algorithms  in  detail  using  one  specific  case  (collective  transfer 
function  at  fj,=0.333),  we  applied  the  algorithms  to  estimate  the  HTFs  of  the  ATR  system  (collective, 
cosine  cyclic,  and  sine  cyclic  at  ^11=0.14,  0.2,  0.267,  0.333).  Figure  28(a)-(c)  show  the  three  HTFs  of  the 
ATR  system  which  we  believe  physically  represent  the  system  (i.e.,  Go(s),  Gi(s),  and  G2(s))  at  ^=0.333, 
for  collective,  cosine  cyclic,  and  sine  cyclic  input,  respectively.  The  results  in  Figure  28  show  that  the 
magnitude  of  Go(s)  (corresponding  to  LTI  transfer  function)  is  much  larger  than  other  two  HTFs,  at  least 
by  20  dB.  This  result  implies  that  we  can  achieve  significant  amount  of  vibration  reduction  using 
continuous-time  controllers  based  on  LTI  plant  assumption.  This  is  important,  in  that  it  indicates  that 
classical  control  laws  should  work  well  for  the  controlling  the  ATR  rotor,  which  should  greatly  simplify 
the  design  of  closed-loop  controllers  for  the  ATR  rotor. 

Figure  29(a)-(c)  show  the  magnitude  and  phase  of  Go(s)  of  the  ATR  system  at  ^=0.14,  0.2,  0.267, 
0.333,  for  collective,  cosine  cyclic,  and  sine  cyclic  input,  respectively.  The  results  in  Figure  29  show  that 
the  dynamics  of  the  ATR  system  behave  as  expected,  in  that  the  collective  transfer  functions  are  hardly 
affected  by  advance  ratio,  while  the  sine  cyclic  transfer  functions  are  proportional  to  advance  ratio,  both 
in  the  low-frequency  range.  Also,  the  results  indicate  that  the  HTFs  of  the  ATR  system  are  not  so 
sensitive  to  advance  ratio,  especially  their  phase,  which  implies  that  a  controller  carefully  designed  with 
fixed  gains  should  work  adequately  for  various  values  of  advance  ratio. 
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C.2  Closed-loop  Experimental  Results 

Based  on  the  results  of  the  system  identification  for  the  ATR  system,  we  designed  and  implemented 
continuous-time  higher  harmonic  control  (HHC)  algorithms  to  reduce  the  vibration  of  helicopters.  We 
chose  continuous-time  controllers,  rather  than  traditional  discrete-time  HHC  algorithms  developed  in  the 
early  1980’s,  because  our  previous  research  results  show  that  continuous-time  controllers  yield  better 
closed-loop  performance  than  similar  discrete-time  controllers  [4],  More  importantly,  the  results  show 
that  the  discrete-time  controller,  which  implicitly  assumes  that  the  rotor  response  is  quasisteady,  can  drive 
the  rotor  unstable. 

Furthermore,  we  designed  controllers  based  on  the  assumption  that  the  rotor  system  is  linear  time 
invariant  (LTI),  considering  that  the  magnitude  of  Go(s)  is  much  larger  than  other  HTFs.  That  is,  we 
assumed  that  the  rotor  system  behaves  as  an  LTI  system  with  a  transfer  function  Go(s),  and  we  designed 
continuous-time  controllers  so  that  the  loop  transfer  function  (i.e.,  Go(s)  times  the  controller  transfer 
function)  has  a  high  magnitude  at  the  target  frequencies,  while  preserving  enough  gain  and  phase 
margins. 

We  tested  our  continuous-time  HHC  systems  using  various  gain  combinations  for  each  control  input 
(i.e.,  collective,  cosine  cyclic,  and  sine  cyclic),  and  flight  conditions  (;U=0.14,  0.2,  0.267,  0.333).  For  most 
cases,  we  could  reduce  the  vibration  level  significantly  at  the  target  frequencies  (IP  and  4P  components 
were  selected  as  target  frequencies  in  this  experiment).  Table  1  summarizes  the  level  of  vibration 
reduction  we  could  achieve  for  various  closed-loop  experiments  conducted  at  Langley.  Figure  30  shows 
the  typical  closed-loop  performance  obtained  using  our  control  systems,  with  ratios  of  open-loop  response 
to  closed-loop  response.  Note  that  negative  and  positive  ratios  indicate  attenuation  and  amplification  of 
the  vibration,  respectively.  Figure  30  indicates  that  more  than  30  dB  of  reduction  can  be  achieved  at  both 
IP  and  4P  components  using  our  continuous-time  HHC  system. 

This  section  briefly  describes  the  results  of  closed-loop  experiments  conducted  at  Langley.  For  more 
information  about  the  closed-loop  results,  refer  to  [5]  (attached  to  this  report). 
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Figure  30:  Typical  open  &  closed-loop  spectra  of  the  hub  normal  shear  vibratory  load. 
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Table  1:  Results  of  vibration  reduction  for  ATR  system  (cited  from  [5]) 


Case  Name 

Open-loop  RMS 
normal  load  (lb) 

Closed-loop  RMS 
normal  load  (lb) 

Performance 
(Reduction  level,  dB) 

Cycl 

22.35 

13.77 

4.2 

Cyc2 

22.45 

15.14 

3.4 

Cyc3 

22.52 

12.06 

5.4 

Cyc4 

22.55 

8.87 

8.1 

Cyc5 

22.61 

7.19 

10.0 

CollCycl 

22.60 

12.21 

5.3 

CollCyc2 

22.50 

8.75 

8.2 

CollCyc3 

22.26 

7.60 

9.3 

CollCyc4 

10.71 

0.03 

51.1 

CollCycS 

10.63 

0.04 

48.5 

CollCycb 

14.06 

1.41 

20.0 

CollCycV 

14.01 

0.36 

31.8 

CollCycS 

21.45 

10.15 

6.5 

CollCyc9 

21.35 

8.72 

7.8 
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